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ABSTRACT 


The  unsteady  compressible  laminar  boundary  layer  on  a  step- 
wise  accelerated  semi-infinite  flat  plate  at  zero  incidence  is  investi¬ 
gated  for  Prandtl  numbers  different  from  unity.  A  semi -similarity 
transformation  is  employed  to  transform  the  pertinent  boundary  layer 
equations  such  that  the  new  independent  variables  become 
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The  equations  are  uncoupled  for  the  case  of  a  zero  pressure 

gradient  and  are  further  linearized  before  a  numerical  solution  is 

obtained.  The  energy  equation  with  Pr  =  0.72  is  solved  in  two  different 

ways.  First,  an  initial  condition  is  obtained  by  assuming  that 

|^-  (n,0)  =  0,  where  g(n,T)  =  p-  ^n,T^  and  second,  by  employing  an 

e 

empirical  steady  state  enthalpy  distribution  as  a  solution  for  the 
unsteady  case.  For  a  range  of  Mach  numbers  5  £  M  £  10  the  transient 
wall  shear  stress,  heat  transfer  and  weak-interaction  induced  pressure 
are  obtained  for  both  the  isothermal  and  insulated  wall. 
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CHAPTER  I 

STATEMENT  OF  THE  PROBLEM  AND  ASSOCIATED  LITERATURE 


1.1  Introduction 

From  boundary  layer  theory,  one  effect  of  the  viscous  boundary 
layer  on  the  remaining  inviscid  flow  is  an  outward  deflection  of  the 
streamlines  near  the  body  such  that  the  effective  shape  of  the  body 
is  changed  due  to  the  reduced  mass  flow  within  this  boundary  layer. 

A  measure  of  this  displacement  of  the  external  flow  is  given  by  the 
displacement  thickness  defined,  for  compressible  flows,  as  [1,2]* 


(i.i) 


For  slender  bodies,  such  as  a  flat  plate,  moving  at  sub¬ 
sonic  or  low  supersonic  speeds,  6*  is  of  the  order  of  1//Re  so  that 

X 

for  high  local  Reynolds  numbers  the  manifestations  of  streamline  dis¬ 
placement,  such  as  induced  pressures  for  example,  can  usually  be  ne¬ 
glected.  For  the  same  body  flying  at  hypersonic  speeds,  however,  the 
temperature  in  the  boundary  layer  will  be  extremely  high  due  to  frictional 
heating  and  consequently  the  density  v/ill  be  very  low.  Thus,  as  equation 

(1.1)  dictates,  the  displacement  thickness  will  be  large  and  in  fact, 

2 

the  resulting  streamline  deflection  will  be  of  the  order  of  M  //Re 


*Numbers  within  square  brackets  refer  to  corresponding  references. 
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and  the  induced  pressure  due  to  the  interaction  of  this  thick  boundary 
layer  with  the  Mach  waves  in  the  free  stream  is  of  the  order  of 
M0  //Rex  [3].  It  becomes  clear,  then,  that  at  hypersonic  speeds  viscous- 
inviscid  interactions  must  be  taken  into  account  and  aside  from  obtain¬ 
ing  the  magnitudes  of  the  induced  pressure  it  is  desirable  to  deter¬ 
mine  what  effect  such  interaction  would  have  upon  the  surface  skin 
friction  and  heat  transfer  coefficients. 

In  order  to  study  and  formulate  the  above  phenomenon  the 
following  hypothetical  problem  is  posed.  Figure  1.1  depicts  a  semi- 
infinite  flat  plate  with  a  "sharp"  leading  edge  moving  in  air  at 
hypersonic  speeds.  The  word  "sharp"  is  used  in  the  same  sense  as  in 


Fig.  1.1  Flat  Plate  at  Hypersonic  Speed 
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Reference  [3]  and  that  is  to  mean  that  the  leading  edge  of  the  plate 
has  essentially  no  effect  on  the  inviscid  pressure  distribution  along 
the  surface.  A  blunted  slender  body,  for  example,  would  give  rise  to 
another  interaction  phenomenon  usually  referred  to  as  vorticity  inter¬ 
action  but  we  wish  to  concern  ourselves  only  with  the  generally  more 
important  type  of  "self  induced"  interaction  commonly  referred  to  as 
a  pressure  interaction.  The  plate  is  assumed  to  be  at  some  arbitrary 
hypersonic  Mach  number  and  zero  angle  of  attack  with  temperature  and 
velocity  profiles  at  their  corresponding  initial  steady  states.  The 
plate's  velocity,  II,  is  then  impulsively  increased  by  a  small  amount 
of  the  order  of  1%. 

Because  of  the  presence  of  the  boundary  layer  the  external 
flow  near  the  leading  edge  must  turn  as  is  illustrated  by  the  sketch 
of  a  typical  streamline.  The  result  is  an  initial  compression  lead¬ 
ing  to  the  formation  of  a  shock  wave  which,  at  high  Mach  numbers  be¬ 
comes  quite  strong.  The  viscous-inviscid  interaction  comes  into 
play  because  the  flow  expansion  following  the  shock  determines  the 
free  stream  pressure  which  affects  the  redistribution  of  the  boundary  • 
layer  which  in  turn  affects  the  free  stream  pressure.  For  the  steady 
states  approximate  relations  governing  the  free  stream  induced  pressure, 
p  ,  have  already  been  derived  and  these  will  be  discussed  briefly  in 
the  following  section.  Of  prime  importance  also,  however,  are  the 
problems  of  maneuverability  and  control  of  hypersonic  vehicles  and 
these  are  entirely  dependent  upon  the  temporal  variations  of  induced 
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pressures,  shearing  stresses  and  heat  transfer.  The  research  in  this 
direction  has  so  far  been  scarce  and  in  particular  all  available  solu¬ 
tions  have  involved  the  common  assumption  that  the  Prandtl  number  is 
equal  to  unity.  The  mathematical  formulation  and  pertinent  assumptions 
used  in  this  work  are  exactly  the  same  as  found  in  Reference  [4]. 

That  is,  the  two-dimensional,  non-steady,  compressible  boundary  layer 
equations  are  transformed  by  the  use  of  the  Dorodni tsyn-Howarth 
variable  and  the  non-steady  stream  function.  A  semi-simi lari ty  trans¬ 
formation  reduces  the  number  of  independent  variables  from  three  to 
two.  The  resulting  momentum  and  energy  equations  are  uncoupled  by 
considering  only  the  so  called  "weak  interaction"  region  which  is  the 
region  far  downstream  from  the  leading  edge  where  the  free  stream 
velocity  Ue  can  be  considered  a  constant  and  thereby  causing  the 
pressure  gradient  to  vanish.  The  momentum  equation  is  linearized 
using  the  method  of  small  perturbations  and  the  resulting  linear 
equation  is  solved  with  the  aid  of  a  digital  computer  employing  a 
recurrence  type  numerical  procedure  which  was  also  used  in  Reference 
[4].  The  solution  to  the  momentum  equation  obtained  in  this  work  is 
slightly  more  accurate  since  discretization  errors  have  been  reduced 
by  using  more  nodal  values  in  the  computational  molecule,  reducing  the 
step  sizes  in  the  space  and  time  directions  and  employing  double  pre¬ 
cision  for  all  numerical  work. 

The  major  concern  of  the  present  analysis  is  to  obtain  the 
solution  to  the  energy  equation  for  any  arbitrary  Prandtl  number  and 
in  particular  for  Pr  =  .72  which  corresponds  more  closely  to  the  value 
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for  air.  The  solutions  are  obtained  numerically  using  a  similar  re¬ 
currence  technique  as  was  used  for  the  momentum  equation.  The  effects 
of  this  change  on  the  induced  pressures,  skin  friction  and  heat  transfer 
are  discussed  and  compared  with  those  results  for  Pr  =  1  given  in 
Reference  [4].  For  the  case  of  Pr  =  1  a  closed  form  solution  is  avail¬ 
able  both  for  the  insulated  wall  and  the  isothermal  wall  [4].  However, 
for  Pr  f  1  an  initial  condition  on  g  e  h/hg  is  required  for  the  nu¬ 
merical  solution  of  the  unsteady  energy  equation  since  the  condition 
at  infinity  cannot  be  satisfied  due  to  the  fact  that  a  solution  to 
the  unsteady  momentum  equation  cannot  be  developed  beyond  values  of 
t  =  l/f£  .  Thus,  it  will  be  demonstrated  that  by  making  the  reasonable 
assumption  3g/ (n ,0)  =  0,  a  plausible  initial  condition  on  g  can  be 
derived. 

In  Section  3.5  an  alternative  candidate  for  a  solution  to 
the  unsteady  energy  equation  will  be  introduced.  This  solution  will 
be  in  the  form  of  an  empirical  enthalpy  distribution  for  Prandtl 
numbers  differing  only  slightly  from  unity.  The  empirical  distribu¬ 
tion  will  also  be  seen  to  reduce  to  the  closed  form  solutions  mentioned 
above  when  Pr  =  1 . 

1.2  Review  of  Associated  Literature 

The  unsteady  problem  of  a  sudden  increase  in  velocity  of  a 
body  already  in  motion  with  established  velocity  and  temperature  pro¬ 
files  in  the  boundary  layer  is  quite  similar  to  the  problem  of  impulsive 
start  from  rest.  In  fact,  it  has  been  shown  (see  Sears  [5]  for  example) 
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that  at  the  time  of  change  and  except  at  the  solid  boundary  (i.e.  the 
plate  surface  in  our  case)  the  new  velocity  distribution  is  obtained 
simply  by  superimposing  the  suddenly  created  velocity  on  the  established 
velocity  profile  existing  in  the  boundary  layer  just  before  the  change. 
Hence,  some  of  the  arguments  pertaining  to  the  transient  velocity  dis¬ 
tribution  in  the  boundary  layer  of  a  body  starting  impulsively  from 
rest  can  equally  well  be  applied  to  the  problem  under  consideration 
here.  The  difference,  however,  being  that  the  secondary  boundary 
layer,  which  at  the  instant  of  the  impulsive  motion  is  in  the  form  of 
a  vortex  sheet  of  zero  thickness  at  the  plate  surface,  grows  in  thick¬ 
ness  and  interacts  in  a  non-linear  way  with  the  original  or  primary 
boundary  layer  already  in  existence  before  the  impulsive  change  [6]. 
Because  of  the  variation  of  both  velocity  and  temperature  normal 
to  the  plate  within  the  primary  boundary  layer,  it  would  not  be  reason¬ 
able  to  suppose  that  the  entire  time  history  of  the  growth  of  the 
secondary  layer  interacting  with  the  primary  layer  and  ultimately 
reaching  a  new  steady  state,  would  be  quite  similar  to  the  complete 
development  of  a  boundary  layer  from  start  of  motion  to  the  final  steady 
state  resulting  from  impulsive  motion  from  rest. 

The  study  of  boundary-layer  formation  on  an  infinite  flat 
plate  after  an  impulsive  start  of  motion  dates  back  as  far  as  Stokes, 
1901  [2].  The  solution  to  this  problem  for  small  times  (sometimes 
referred  to  as  the  "Rayleigh  problem")  is 

=  U  erf  (— ^-) 
c  2/vt 


u 


(1.2) 


' 


, 
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with 


u (0 ,t)  =  U  . 


The  solution  (1.2)  is  derived  from  the  fact  that  at  the  start  of 
boundary  layer  development  the  effects  of  diffusion  far  outweigh  those 
of  convection.  At  larger  times  from  the  start  of  motion  the  rate  of 
diffusion  balances  the  rate  of  convection  and  the  flow  field  closely 
approximates  a  Blasius  type  of  flow  which  is  independent  of  time. 

J.T.  Stuart  (see  Chapter  VII  of  Reference  [6])  offers  an  enlightening 
analysis  of  bodies  starting  from  rest  from  a  dimensional  analysis 
point  of  view.  For  a  semi-infinite  flat  plate  impulsively  started 
from  rest  it  turns  out  that  t  =  Ut/x  becomes  a  key  variable  in  the 
mathematical  formulation  of  the  flow.  If  t  is  the  elapsed  time  from 
start,  x  the  distance  measured  from  the  leading  edge  of  the  plate  and 
U  the  constant  "step  input"  velocity  then  for  x  very  large  the  flow 
becomes  independent  of  time  and  closely  assumes  the  form  of  a  Blasius 
flow.  This  flow  occurs  for  a  certain  distance  downstream  from,  but 
not  including,  the  leading  edge.  For  x  small  the  flow  is  of  the  Ray¬ 
leigh  type  and  is  not  affected  by  the  leading  edge.  This  flow  occurs 
in  a  region  extending  from  some  point  downstream  from  the  leading  edge 
to  infinity.  A  transition  region  must  therefore  exist  between  the 
two  types  of  dominant  flows  and  this  region  depends  on  x.  As  time 
increases,  the  Blasius  region  extends  farther  and  farther  down  the 
plate  and  as  t  -*  «>  ,  this  region  dominates  everywhere  whereas  for  t 


->  00 


% 


■ 
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very  small,  the  Rayleigh  flow  was  dominant  almost  everywhere. 

Stewartson  [7]  obtained  extreme  solutions  to  the  momentum 
equation  simplified  to 


3u  ,  ,,  3u 
Ft  +  Ue  3x 


=  V 


32u 

sy2 


(1.3) 


where  Rayleigh's  analogy  was  used  to  replace  the  convection  terms 


u~-+v— by  U 
3x  3y  J  e  8x 


The  solutions  obtained  are 


u  =  U  erf  — }  t  <  1 
e  2A)t 


and 


u  =  Ue  erf  {2  /  w}  X  >  1 


(1.4) 


Strictly  speaking  these  solutions  are  valid  only  near  the  boundary 

layer  edge  where  u  *  U  but  they  do  show  qualitatively  the  motion  of 

e 

the  Blasius  to  Rayleigh  transition  region  down  the  plate  as  time  in¬ 
creases.  Using  the  momentum- integral  method  as  an  alternative  approxi¬ 
mation,  Stewartson  obtained  for  the  shear  stress  at  the  wall 

(yfj-)  =  0.534  pU 
y=0 


x  <  2.65 


(1.5) 


' 
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(u  f;)  =  0.328  p  U  /  ^  t>2.65.  (1.6) 

y  y= 0  e  A 

Again  an  abrupt  change  from  a  time  dependence  to  an  x-dependence  at 
some  definite  value  of  t  is  indicated  which  does  not  seem  physically 
reasonable.  It  v/as  further  shown  by  Stewartson  that  a  power  series  in 
t  to  account  for  gradual  decreasing  time  dependence  and  increasing 
x-dependence  cannot  be  constructed  and  this  is  also  indicated  in  the 
discussion  by  Stuart. 

Many  investigators  including  S.I.  Cheng,  F.K.  Moore,  Simon 
Ostrach  and  others  have  studied  similar  problems  and  a  resume  of  their 
works  may  be  found  in  Reference  [4].  In  each  case  a  singularity  was 
encountered  in  the  neighbourhood  of  t  =  1.  Rodkiewicz  and  Reshotko  [4] 
solved  numerically  the  unsteady  momentum  equation  describing  the  present 
problem  but  the  solution  could  not  be  generated  very  far  beyond  x  =  1. 
Rodkiewicz  and  Gupta  [8]  extended  the  solution  slightly  to  x  =  1.3 
by  use  of  double  precision  but  the  apparent  success  is  misleading 
since  the  present  writer  has  solved  the  same  equation  using  double 
precision  as  well  but  reducing  discretization  errors  by  including  extra 
nodal  values  in  the  computational  molecule  and  reducing  both  stepsizes 
in  space  and  time  to  1/1 0th  of  the  values  used  in  both  [4]  and  [8]. 

The  results  show  a  discontinuity  beginning  at  a  value  of  x  only  slightly 
greater  than  unity.  However,  no  serious  attempt  has  been  made  in  this 
work  to  resolve  this  problem. 

The  effects  of  Prandtl  number  variation  on  steady  laminar 


«r 


- 

- 
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compressible  flows  have  been  quite  thoroughly  discussed  by  several 
authors.  Fliigge-Lotz  and  Arlo  F.  Johnson  [9],  for  example,  concluded 
that  the  drag  of  an  insulated  flat  plate  at  zero  angle  of  attack  is 
unaffected  by  Mach  and  Prandtl  numbers  and  that  Prandtl  number  vari¬ 
ation  appears  to  create  a  significant  change  only  in  the  enthalpy 
profile.  They  do  not,  however,  consider  boundary-layer  shock-wave 
interaction.  For  unsteady  flows  the  results  of  Ostrach  [10],  for 
example,  are  for  a  Prandtl  number  of  .72  but  do  not  lend  themselves 
for  direct  comparison  with  this  work  since  the  plate  is  assumed  to 
have  a  continuous  time-dependent  velocity.  For  the  problem  to  be 
considered  in  this  research,  the  author  has  not  been  able  to  find  any 
solutions  for  Prandtl  numbers  different  from  unity. 

Estimates  of  the  steady  state  induced  pressure  are  usually 
obtained  with  quite  accurate  results  from  the  "tangent-wedge"  method 
[3,11].  The  resulting  expression  is 


1  +  y  M 

00 


d6j 

dx 


y(y+i) 

4 


(M 

'  oo 


d6*< 
dx  ' 


(1.7) 


with  M  >  1  and  M  (dS*/dx)  «  1. 

00  00 

Thus  for  the  weak  interaction  case  one  might  simply  say  Pg/p^  ~  1  as 
a  first  approximation.  Such  a  simplifying  approximation  would,  how¬ 
ever,  defeat  the  basic  purpose  of  an  interaction  study  and  hence  we 
will  concern  ourselves  with  the  transient  contribution  to  the  induced 
pressure.  In  Moore  [12]  the  outward  displacement  of  the  boundary 
layer  is  likened  to  a  mild  constant-pressure  explosion  after  the  in- 


. 
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stant  of  the  impulsive  change  and  associated  with  the  resulting  "sound" 
waves  sent  out  into  the  inviscid  stream  is  a  pressure  rise  which  from 
the  acoustic  approximation  is  given  by  p  a  v  . 

c  0  0 

The  acoustic  approximation  may  be  derived  from  Lighthill's 
piston  theory  [13]  as  discussed  in  a  paper  by  J.W.  Miles  [14].  The 
perturbation  pressure  may  be  regarded  as  resulting  from  the  motion  of 
a  piston  normal  to  the  plate  and  this  disturbance  may  be  treated  as 
a  simple  wave.  The  perturbation  pressure  is  then  given  by  the  relation 

Apv  =  pav1  (1.8) 


where  p  and  a  are  the  local  density  and  speed  of  sound  and  v'  is  the 
unsteady  perturbation  velocity  normal  to  the  perturbation  wave  front. 
The  pressure,  density  and  sonic  speed  at  the  piston  thus  depend  only 
on  the  instantaneous  piston  velocity,  W  say,  and  are  given  by 

r 


2 


_L  = 


(^)Y 


-1 


JL 

Pe 


=  (~) 


a_xy-l 


a_ 

a_ 


i  W 

1  +  ’r  ^ 


(1.9) 


Writinq  W  =  W  +  v 1  wi th  W  =  U  0  and  0  representing  the  local  steady 
3  p  e 

inclination  of  the  surface  to  the  free  stream,  and  using  W  =  W  in 
(1.9)  for  the  calculation  of  p  and  a,  we  obtain 


apv  «  Pe  ae  V  [l  +  ^  Me  eF1 


t 
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=  Pe  ae  v'0  +  +  Me2  ®Z  +  °(Me3e3)]  •  O-TO) 

Since  for  our  problem  0  =  d6*/dx  <<1  and  v'  =  v  expression  (4.11) 

e 

for  the  induced  pressure  follows.  The  normal  component  of  velocity 
v  is  obtained  from  the  expression  for  the  unsteady  displacement 
thickness  given  by  Moore  and  Ostrach  in  Reference  [15],  and  the  re¬ 
sult  is  an  expression  for  the  transient  contribution  to  the  induced 
pressure  given  by  equation  (4.28). 

The  usual  Prandtl  boundary  layer  assumptions  have  been  adopted 
in  the  works  cited  above  and  will  also  be  used  in  the  present  research. 
Although  the  boundary  layer  in  hypersonic  flow  will  be  quite  thick  it 
has  been  pointed  out  by  Shen  [16]  that  the  assumptions  remain  valid 
for  6*/x  sufficiently  small.  For  the  weak  interaction  region,  x  will 
be  large  so  that  this  condition  will  be  met  in  the  present  work.  Also, 
it  is  shown  by  Butler  [17]  from  a  detailed  investigation  that  the 
pressure  across  the  hypersonic  boundary  layer  tends  to  vary  less  (i.e. 
remain  fairly  constant)  at  increasing  values  of  x,  where  x  is  measured 
from  the  leading  edge  of  the  plate. 

Finally  we  note  that  the  boundary  layer  under  investigation 
is  assumed  to  be  laminar.  For  a  typical  cruise  altitude  of  a  hypersonic 
vehicle,  say  90,000  ft.  slightly  above  the  stratosphere,  the  free 
stream  temperature  is  approximately  -  58°F  (402°r).  The  corresponding 
sonic  speed  is  then  approximately  985  ft/sec.  Hence  the  free  stream 
velocity  for  a  Mach  number  of  10  would  be  9850  ft/sec.  This  would  give 
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a  local  Reynolds  number  based  on  a  nominal  distance  of  10  ft.  from  the 

g 

leading  edge  of  about  9.85  x  10  which  is  by  far  above  the  critical 
Reynolds  number  for  a  flat  plate.  Other  effects  on  the  laminar  to 
turbulent  transition  such  as  pressure  gradient,  for  example,  must 
also  be  considered.  In  our  case  the  effect  of  heat  transfer  and 
compressibility  is  a  favourable  one  since  it  is  known  that  transfer 
of  heat  from  the  boundary  layer  to  the  wall  exerts  a  stabilizing 
influence  by  increasing  the  critical  Reynolds  number  [2].  In  fact, 
there  appear  to  be  certain  combinations  of  Mach  number  and  wall  en¬ 
thalpies  (g  )  for  which  the  boundary  layer  is  completely  stable  i.e. 
w 

for  which  the  critical  Reynolds  number  becomes  infinite.  We  therefore 
retain  our  mathematical  model  on  the  assumption  of  a  laminar,  com¬ 
pressible,  unsteady  two-dimensional  boundary  layer. 


•*  ' 
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CHAPTER  II 


THE  GOVERNING  BOUNDARY  LAYER  EQUATIONS 

2.1  Unsteady  Prandtl  Boundary  Layer  Equations 

The  general  form  of  the  unsteady,  compressible,  laminar, 
two-dimensional  boundary  layer  equations  with  the  adoption  of  Prandtl's 
boundary  layer  assumptions  for  hypersonic  flow  are  as  follows  [12,18]: 


(2.1) 


(2.2) 


(2.3) 


The  boundary  conditions  are: 


at  y  =  0 


u  =  v  =  0 


h  =  hw(x,t) 


(2.5) 


at  y  =  ye 


u  =  Ue(x,t) 
h  =  he(x,t) 


I 
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The  pressure,  being  assumed  independent  of  y,  imposed  on  the  plate  is 
given  by  its  free  stream  value  according  to  the  solution  of 


1 


3P, 


Pe 


3U  3U 

— 6  +  y  — - 

3t  e  3x 


(2.6a) 


which  is  the  Eulerian  equation  of  motion  [19]  and  follows  immediately 
from  (2.2)  when  specialized  to  the  inviscid  free  stream  flow.  Similarly, 
specializing  (2.4)  to  the  free  stream  one  obtains 


P  J 


3h 


ev3t 


„  3h  3p  3p 

+  U.  ~  +  u  ‘ 


e  3x 


3t 


e  3x 


(2.6b) 


where  U  (x,t)  represents  the  prescribed  non-steady  potential  motion, 
e 

2.2  The  Dorodni tsyn-Howarth  Variable  and  the  Non-Steady  Stream  Function 
It  is  possible  to  rewrite  the  equations  (2.1)  through  (2.4) 
in  a  form  similar  to  that  for  the  incompressible  boundary  layer  by 
adopting  the  Dorodni tsyn-Howarth  variable  [18] 

y 

Y  =  dy  (2.7) 

J  pref 
o 

and  introducing  the  stream  function  ip(x,y,t)  so  that  the  velocity 
components  are  related  to  it  by  the  equations 

pref  3^  _  3£_ 
p  3y  3Y 


u 


(2.8) 
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v 


Pref  / .  3Y n 
p  '  3x  3t; 


(2.9) 


The  equation  of  continuity  (2.1)  is  then  automatically  satisfied. 

Because  of  the  high  temperature  rise  in  the  boundary  layer 
at  high  Mach  number  flows,  it  becomes  necessary  to  take  into  account 
the  effect  of  temperature  on  the  properties  of  the  gas  and  in  particu¬ 
lar  its  viscosity.  An  accurate  relation  between  viscosity  and  temper¬ 
ature  is  given  by  Sutherland's  equation 


y 

^ref 


=  G 


T  ^  Tref  +  S 


ref 


T  +  S 


(2.10) 


where  S  is  Sutherland's  constant  which  for  air  is  approximately  110°K  [2]. 
Following  Chapman  and  Rubesin  [20]  equation  (2.10)  may  be  approximated 
by  the  relation 


—it—  =  C(x)  =5-3—  (2.11) 

^ref  ref 


where  C(x)  is  fixed  so  that  the  viscosity  is  accurately  determined  in 
the  more  important  region  near  the  surface  of  the  plate.  If,  for 
example,  equations  (2.10)  and  (2.11)  are  matched  at  a  temperature  cor¬ 
responding  to  the  average  surface  temperature,  I w  ,  then  [20] 


+  S 
+~S“ 


C(x)  =  (/ 


(2.12) 


. 


Finally,  assuming  a  perfect  gas  with  thermal  equation  of  state 
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p  =  pRT 


(2.13) 


the  viscosity  relation  becomes  as  in  [12] 


irVc{x>A  .  (2.14) 

Mref  H 

Substituting  equations  (2.7),  (2.8),  (2.9)  and  (2.14)  into  equations 
(2.2)  and  (2.4),  the  momentum  and  energy  equations  become  respectively 

[4] 


d2\p  3i P  92i p  9i|j  92 p  _  pe  ..  9Ue  r  a3i|) 

3Y3 1  9Y  9x9Y  "  9x  ^ 2  p  e  9x  L  vref  y3 


(2.15) 


ah.  ,  a^  ah  9^  9h_  = 

9t  T  9Y  9x  “  9x  9Y 


9U  r  A  A 

p  9Y  Me  3>r  +  vref  Pr  +  vref  (2-16) 


e  dip 


U 


It  has  been  further  assumed  that  the  Prandtl  number  Pr  is  constant  for 


air, 


2.3  Boundary  Layer  Equations  in  Similarity  Form 

As  demonstrated  by  Rodkiewicz  and  Reshotko  [4]  equations 
(2.15)  and  (2.16)  may  be  further  simplified  by  the  introduction  of  two 
similarity  variables  which  reduce  the  number  of  independent  variables 
from  three  (x,Y,t)  to  two  (ti,t).  The  following  transformation  of 


variables  was  used: 
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n 


(m+1)  Ue 
Fvref  Cx 


)Y 


T 


u  t 


(2.17) 


(2.18) 


where  the  free  stream  velocity  must  be  in  the  form  of  a  power  law 


Ue  =  exm  .  (2.19) 

A  dimensionless  stream  function,  f(n,x),  becomes  the  new  dependent 
variable  in  the  momentum  equation  and  is  defined  by 

2^rpf  ^  U 

'  »  =</  -  %■!)  e-)f(n,T)  .  (2.20) 

The  energy  equation  is  transformed  by  letting 


h  =  g(n,T)  G(x,t) 


(2.21) 


where  g(n,T)  becomes  the  new  dimensionless  enthalpy  function.  Upon 
substitution  of  the  above  expressions  into  equations  (2.15)  and  (2.16) 
one  obtains 


.  U  9  .  _  .  BU  , 

f'  f  +  f  IT  -  Uef,f’  7  +  f'f'T  33T  ‘  27  ff'  'Ue  +  fUef' '  7 


-  f  f  '  '  T 


i  ffl,  ^  jl!^ 
8x  2 


3x  -  h  3X  '  CVffl2Tf"  "  0 
e 


(2.22) 
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and 


gG^+gfl^G^f  +xmf  Ggf  +g|fUef 


,  U  ,  U  U  .  u 

{  -f  Gfg'  -  j  m  -f  fGg '  +  t  -f  f Gg '  -  tfGg'm  -j 


U 


U 


f  UD2f'm  4  +  Gg' 1  +  -f-  f  '2 


(2.23) 


where 


m+1 


=  8f 
“  3t 


f 


=  if 
3n 


It  is  assumed  [4]  that  at  hypersonic  speeds  the  free  stream 
adjusts  itself  instantaneously  to  the  new  state  resulting  after  any 
sudden  change  in  the  free  stream.  Mathematically  this  means  that 


3h  3p 
_ e  _  re 

8t  "  3t 


(2.24) 


Consequently  expressions  (2.6a)  and  (2.6b)  can  be  combined  to  give 


e  3x 


(2.25) 
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and  using  (2.19) 


and  letting  G(x,t)  =  G(x)  =  h  (x)  one  obtains  finally 

0 


3G 

3x 


(2.26) 


(2.27) 


Substitution  of  the  pertinent  expressions  into  equations  (2.22)  and 
(2.23)  yields  the  equations  in  their  final  similarity  form,  namely 


(2-g)f 1  +  2(1-B)t{ f f ' ' -f 1  f ' )  -  f'"  -  ff"  -  g(g-f'2)  =  0  (2.28) 


and 


Pr(2-g)  [1  -  2  g  -  g"  -  Pr[f  -  2(l-g)Tf]g‘ 


-  Pr(y-l)  M02f' ,2  =  0 


(2.29) 


where 


=  (y-1)  Me2  . 


and 


- 
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For  the  case  under  consideration  in  this  work  (zero  pressure 
gradient,  m  =  3  =  0)  equations  (2.28)  and  (2.29)  become  uncoupled  and 
reduce  to 


2f '  +  2x(ff 1  1  -f 1 f 1 )  -  f1"  -  ff"  =  0  (2.30) 

2Pr(l -xf 1 )  g  -  g"  -  Pr(f-2xf)  g'  -  Pr(y-l)  Me2f"2  =  0  (2.31) 

For  the  final  steady  state,  equation  (2.30)  becomes  the  familiar 
Blasius  equation 

f2 '  '  +  V2*  =  0  (2.32) 


and  equation  (2.31)  becomes 

92 ‘  +  Pff232  +  My-1)  M02f“2  =  0  .  (2.33) 

2.4  Linearization  of  the  Momentum  Equation 

Equation  (2.30)  is  a  non-linear  partial  differential  equation. 
The  method  of  small  perturbations ,  often  used  to  linearize  such 
equations,  is  in  general  not  applicable  to  equations  governing  the 
flows  in  the  hypersonic  regime.  Since  we  are  considering  only  a  very 
small  velocity-step  increase  of  the  plate,  i.e.  a  1%  velocity  jump 
for  only  low  hypersonic  Mach  numbers  (maximum  =  10),  it  has  been 
decided  to  adopt  this  linearization  technique  as  was  done  in  Reference 
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[4].  Furthermore,  Gupta  and  Rodkiewicz  [21]  have  solved  the  non-linear 
momentum  equation  (2.28)  for  the  strong  interaction  region  using  a 
numerical  procedure  developed  by  Smith  and  Clutter  [2].  Their  results 
show  that  the  maximum  discrepancy  in  the  shear  stress  parameter  f ' 1 
from  a  linearized  solution  to  (2.28)  is  less  than  1%  for  a  1%  jump 
in  the  free  stream  velocity.  We  therefore  assume  that 

f(n,x)  =  f2(n)  +  Af(n»T)  (2.34) 

where  f2(n)  is  the  known  solution  to  the  Blasius  equation  (2.32). 

(See  for  example  Reference  [6]). 

From  (2.34)  we  have 


f'  =  fxz  +  (Af)' 
f"  =  f“  +  (Af)" 
f'"  =  f2"  +  (Af)'" 

9  • 

f  =  (Af) 

f'  =  (Af)' 


and  substitution  into  (2.30)  yields  the  linear  equation 
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2(Af)'[l-Tf£]  +  2xf‘,(Af)  -  ( Af ) 1  1 1  -  f2(Af)"  -  f“(Af)  =  0  (2.35) 


where  all  higher  order  terms  in  Af  have  been  neglected. 

2.5  Associated  Boundary  and  Initial  Conditions 

When  a  body  is  impulsively  set  into  motion  from  rest  in  a 
viscous  fluid  the  boundary  layer,  immediately  after  the  start  of 
motion,  is  in  the  form  of  an  infinitely  thin  vortex  sheet  which  is  a 
consequence  of  the  non-slip  condition  on  the  surface  of  the  plate. 

In  the  initial  stages  of  growth  of  this  boundary  layer,  downstream 
from  the  leading  edge,  the  effects  of  convection  are  negligible  com¬ 
pared  with  those  of  diffusion  so  that  in  the  mathematical  formulation 
of  the  flow,  equation  (2.2)  may  be  simplified  to 


P 


3u  3p  ,  3  /  3Ux 

n  =  '  a!  +  W  {v  W 


(2.36) 


which  is  linear  and  indicates  a  balance  of  pressure,  inertia  and 
viscous  forces. 

The  physical  problem  we  are  considering  in  this  work  is 
closely  related  to  boundary  -  layer  flow  from  rest.  The  free  stream 
velocity  Ue  is  suddenly  changed  from  one  value  UeQ  to  another  by 
a  small  increment  of  VI  defined  by 


.01 


\ 


(2.37) 


»  ■ 


, 
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Sears  [5]  has  shown  that  in  this  case  the  initial  motion,  since  the 
equation  governing  the  initial  motion  is  linear,  can  be  represented 
by  the  superposition  of  the  initial  flow  pattern  described  above  on 
the  established  flow  pattern  existing  just  before  the  change.  That 
is,  except  at  the  plate  surface,  the  initial  boundary  layer  reacts  as 
if  it  were  inviscid  since  the  sudden  change  has  not  allowed  sufficient 
time  for  viscosity  to  counteract  the  increase  in  velocity  so  that  the 
governing  equation  in  the  major  portion  of  the  original  boundary  layer 
is 


3u  1  3p 

3t  "  p  3x 


(2.38) 


(See  also  Reference  [6]).  The  situation  may  perhaps  be  easier  visua¬ 
lized  with  reference  to  Fig.  2.1. 


Fig.  2.1  Laminar  Velocity  Profiles 
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The  associated  boundary  and  initial  conditions  pertinent  to 
equations  (2.35)  and  (2.31)  are  therefore  as  follows: 

When  y  =  ye,  Y  =  Ye>  n  =  ne,  u  =  Ue2,  f '  (ne ,t)  =  f^Cng.x)  =  1, 

and  Af 1 (n  ,t)  =  0  (2.39) 

When  y  =  0,  Y  =  0,  n  =  0,  u  =  0,  f'(0,x)  =  f^Ojx)  =  0  , 

and  Af 1 (0,x)  =  0  (2.40) 

When  y  =  0,  Y  =  0,  n  =  0,  4)  =  0,  f(0,x)  =  ^(O.t)  =  0 

and  Af(0,x)  =  0  (2.41) 

When  t  =  0,  x  =  0,  u  =  uQ  +  Uel  -  Ug0  and  for  n  >  0 

f 1  (t,>0 ,0)  =  f'  +  e(l-  f'Q) 

so  that 

When  y  =  ye,  Y  =  Ye,  n  = 


Af'(n>0,0)  =  e(l-f^) 


(2.42) 
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and  g(ne,x)  =  1  (2.43) 

When  y=0,Y  =  0,r)  =  0,h  =  hw  and 

g(O.x)  =  gw  (2.44) 

whereas  for  the  case  of  no  heat  transfer  at  the  wall 

g'(0,T)  =  0  (2.45) 


When  t  =  0,  x  =  0,  h  =  h-j  and 

hi 

g(n,0)  =  tJ-  .  (2.46) 

ne2 

An  observation  made  by  Rodkiewicz  and  Reshotko  [4]  concerns 
the  relationship  between  variables  referenced  to  the  initial  and 
final  steady  states.  Recalling  expression  (2.17)  we  have 

(m+l)  U  0 

Tl0=</  2v~~VCx  )  Y  [2. Ala) 

ref 

and 

(m+1 )  U  ? 


(2.47b) 


% 
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so  that 


^0  ,HeO 

n2  ue2  • 

Since  fQ  is  a  function  of  nQ  and  f2  a  function  of  n2  »  boundary  condi¬ 
tion  (2.42)  should  strictly  speaking  be  written  as 

Af' (n2>0,0)  =  e[l-fy/Mn2)]  +  f^(/ JeO  n2)  -  f£(n2)  • 

e  2  e  2 

Considering  equation  (2.37),  however,  we  have 

U  n 

/  ~  =  /(l  -  e)  =  .995 
Ue2 

so  that  rig  =  -995  n2  f°r  a  ^  velocity  jump.  This  small  discrepancy 
would  appear  to  be  insignificant  in  the  numerical  solution  of  the 
equations  and  has  therefore  been  ignored.  Even  for  a  10%  jump  where 
then  rig  =  .95  n2  >  it  would  seem  that  the  difference  is  not  likely  to 
affect  the  results  a  great  deal. 


1 


' 
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CHAPTER  III 

SOLUTIONS  TO  THE  GOVERNING  EQUATIONS 

3.1  Numerical  Solution  of  the  Momentum  Equation 

The  linearized  momentum  equation  (2.35)  subject  to  the  bound¬ 
ary  and  initial  conditions  described  in  Section  2.5  has  been  solved  by 
Rodkiewicz  and  Reshotko  [4]  using  the  technique  of  finite  differences 
and  a  recurrence  relationship  which  is  essentially  a  modification  to 
the  method  of  Gaussian  elimination  [22].  The  range  of  integration  of 
(2.35)  in  the  direction  of  t  extends  from  zero  to  infinity,  however, 
it  turns  out  that  the  numerical  solution  could  be  generated  for  values 
of  t  ranging  from  zero  to  only  slightly  greater  than  1.  The  discon¬ 
tinuity,  although  abrupt  mathematically,  is  a  manifestation  of  a  gradual 
transition  from  a  Rayleigh  type  of  flow  to  a  Blasius  flow  [6].  Tokuda 
[23]  has  considered  a  suggestion  made  by  Van  Dyke  who  postulated  that 
the  mathematical  singularity  might  be  avoidable  by  a  simple  change  of 
coordinate  system  such  as  Euler's  transformation  given  by 


Incorporating  (3.1)  in  equation  (2.35)  yields  the  transformed  momentum 
equation 


I  .  i 


- 


29 


-  Af'"  -  f2Af"  -  f“Af  =  0  . 


(3.2) 


Unfortunately  the  author  has  found  that  the  numerical  solu¬ 


tion  of  (3.2)  could  also  not  be  extended  very  far  beyond  a  value  of 
x  =  1(C  =  1/2).  Now  since  the  range  of  integration  of  the  differential 
equations  in  the  x  direction  extends  from  zero  to  infinity,  it  was  de¬ 
cided  to  use  the  Euler  transformation  (3.1)  for  the  purpose  of  econo¬ 
mizing  in  computer  time  and  storage  since  the  corresponding  range  for 
^  is  0  to  1.  Furthermore,  an  improved  accuracy  in  the  numerical  re¬ 
sults  could  be  realized  since  it  would  be  possible  to  reduce  the  step- 
size  T  in  the  time  direction.  The  energy  equation  (2.31)  in  terms 
of  the  new  time  variable  £  then  becomes 


-  (y-1)  Me2f"2  =  0 


(3.3) 


where  the  dot  over  any  variable  is  now  interpreted  as  differentiation 
with  respect  to  £• 

The  computational  molecule  used  in  conjunction  with  a  finite 
difference  representation  of  derivatives  is  shown  schematically  in 
Fig.  3.1.  It  was  employed  both  in  Reference  [4]  and  [8]  and  utilizes 


. 


. 
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seven  nodal  values  in  a  rectangular  grid. 


-6 — 6- 

T - Y 

-0 - O- 


m  - 1  m 


n  +  | 

n 


n-  I 


/7-  2 


Fig.  3.1  7-Point  Computational  Molecule 

In  conjunction  with  this  molecule  the  finite  difference  representation 
of  the  derivatives  of  equation  (3.2)  are  as  follows: 


||^  =  j  (Af (m,n)  -  Af(m-1 ,n)}  +  0(T)  (3.4) 

|^  =  i  {Af(m,n+1)  -  Af(m,n-1)}  +  0(H2)  (3.5) 

2 

{Af{m,n+1)  -  Af(m,n-1)  -  Af(m-l,n+l) 

+  Af(m-1 ,n-l)}  +  0(T)  +  0(H2)  (3.6) 


\  {Af(m,n+1 )  -  2Af(m,n)  +  Af(m,n-1)}  +  0(H2)  (3.7) 

an  IT 
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{Af(m,n+1)  -  3Af(m,n)  +  3Af(m,n-l) 

an 


-  Af(m,n-2)}  +  0(H)  .  (3.8) 


It  should  be  mentioned  that  a  Taylor  series  expansion  has  been  employed 
in  order  to  eliminate  the  nodal  point  (m,n+2)  for  the  difference  re¬ 
presentation  of  the  third  order  derivative,  equation  (3.8).  If  now 
the  expressions  (3.4)  to  (3.8)  are  substituted  into  (3.2)  the  resulting 
difference  equation  is  of  the  form 


a  s  , .  +  b  s  +  c  s  i+d  s  0  =  e  „  (3.9) 
m,n  m,n+l  m,n  m,n  m,n  m,n-l  m,n  m,n-2  m,n 


where  sm  n  =  Af(m,n)  and 


a  =  1  +  HTf2(m’n)  '  (3.10a) 


b  „  =  -  3T  -  2HTf,(m,n)  +  H3TfV(m,n)  -  2£  ( 1 -5m)  H3f''(m,n)  (3.10b) 
m,n  c.  c.  m  m  c 


cm  n  =  3T  +  HTf2(m,n)  +  H2(l-Cj  [l-Sjl+f^m.n) )]  (3.10c) 


dm  „  =  -  T 
m,n 


(3.1 Od) 


em  n  =  H2 (1  -5m) Cl -Cm(  1  +  (m . n ))  1111s (m-l  .n-1 )  -  s(m-l,n+l)] 


-  2H3  5m(l-5m) (m.n)  s(m-l.n) 


(3.1 Oe) 
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The  numerical  solution  consists  of  a  step  by  step  procedure 
of  finding  's'  completely  along  one  time  line  before  proceeding  to  the 
next  beginning  with  the  known  (or  assumed)  initial  conditions.  That 
is,  the  finite  difference  representation  in  the  time  direction  is  a 
backward  difference  as  is  evident  in  (3.4).  The  nodal  values  on  each 
time  line  are  evaluated  from  the  recurrence  relationship 


s  =  R  s  , +  C 
m,n  m,n  m,n+l  m,n 


(3.11) 


which  follows  from  expression  (3.9)  with 


-  a 


R 


m,n 


m,n  Tb  +  c  R  +  d  R  T~R  ZJ 
v  m,n  m,n  m,n-l  m,n  m,n-l  m,n-2' 


(3.12) 


and 


m,n 


e  -  ( c  C  n+dR  0C  ,+dC  0) 
m,n  v  m,n  m,n-l  m,n  m,n-2  m,n-l  m,n  m,n-27 


Tb  +  c  R  -i  +  d  R  ,R  0i 

v  m,n  m,n  m,n-l  m,n  m,n-l  m,n-2i 


(3.13) 


Beginning  with  the  known  boundary  conditions  at  the  outer  edge  of  the 
boundary  layer  and  at  the  plate  surface,  the  values  of  Rm  n  and  n 
are  calculated  by  marching  up  from  the  plate  surface  to  the  free  stream. 
From  the  boundary  conditions  at  the  wall  and  using  a  three  point  forward 
difference  approximation  to  9s/3ri  it  can  be  shown  that 


1 

4 


(3.14) 


' 


- 


33 


At  the  free  stream  where  n  =  N  say,  the  values  of  s  are  obtained 

m,n 

recursively  marching  downward  using  (3.11)  with  the  starting  value  ob¬ 
tained  from  the  boundary  condition  at  the  outer  edge  of  the  boundary 
layer,  namely 


Tn,N-2 


m,N-l  4-R 


m,N-2 


-  R_ 


4-Rm,N-2  m’H-] 


(3.15) 


It  will  be  noticed  that  the  backward  difference  approxima¬ 
tions  (3.4)  and  (3.6)  have  a  truncation  error  in  the  time  direction 

of  order  T  whereas  the  error  in  the  other  approximations  are  of  order 

2  -2 

H  .  Since  the  maximum  numerical  value  of  Af  is  only  of  order  10  it 

is  imperative  that  the  error  in  the  numerical  solution  be  made  as  small 
as  possible.  Rodkiewicz  and  Reshotko  [4]  have  tried  a  symmetrical 
11-point  computational  molecule  but  with  little  success.  The  author, 
however,  has  successfully  used  a  10-point  molecule  shown  in  Fig.  3.2. 


n 


-0 — 6 — 6-  n  *  \ 


-6 — 6 — 6-  n 


O - 6 - 6-  n-  I 


y-  n-  2 


m-\  m  /77+ 1 


Fig.  3.2  10-Point  Computational  Molecule 
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The  procedure  employed  was  that  only  the  first  time  line  was 
calculated  using  expressions  (3.4)  to  (3.8)  relating  to  the  7-point 
molecule  in  Fig.  3.1.  All  subsequent  time  lines  were  then  evaluated 
using  the  10-point  molecule  of  Fig.  3.2  with  the  following  finite 
difference  approximations. 

f| -  =  27  {3Af(m+l ,n)  -  4Af(msn)  +  Af(m-l,n)}  +  0(T2)  (3.16) 

|7=  {Af(m+l,n+l)  -  Af(m+1  ,n-l )}  +  0(H2)  (3.17) 

2 

Infl’ =  4HT  {3Af('n+1  >n+1 )  "  4Af(m,n+l)  +  Af(m-l,n+l) 

-  3Af (m+1  ,n-l )  +  4Af(m,n-l)  -  Af(m-l,n-l)}  +  0(T2)  +  0(H2)  (3.18) 

2 

\  (Af (m+1  ,n+l )  -  2Af(m+l,n)  +  Af(m+l,n-l)}  +  0 ( H2)  (3.19) 

an  hr 


{Af (m+1  ,n+l )  -  3Af  (m+1  ,n)  +  3Af(m+l  ,n-l ) 

3ri 

-  Af (m+1 ,n-2) }  +  0(H)  (3.20) 

A  three-point  Newton  backward-difference  [2,24]  has  been  used  in  (3.16) 
and  (3.18)  to  approximate  the  derivative  with  respect  to  £  resulting 
in  a  balanced  system  for  a  square  grid  with  truncation  errors  of  0 (T  ) 


' 
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2 

and  0(H  )  in  the  time  and  space  directions  respectively.  With  appro¬ 
priate  revisions  to  quantities  (3.10a)  to  ( 3 . 1 Oe )  the  numerical  solu¬ 
tion  proceeds  in  an  analogous  manner  as  described  previously  except 
that  now  the  recurrence  relation  (3.11)  is  applied  to  each  (m+l)st 
time  line  in  succession.  A  further  refinement,  described  in  Section 
3.2,  was  introduced  before  the  numerical  solution  was  carried  out  and 
this  involved  the  reduction  of  the  step  sizes  H  and  T  to  one-tenth  of 
the  size  used  in  References  [4]  and  [8].  Furthermore,  all  calculations 
were  done  in  double  precision  on  the  IBM  360.  The  solution  curves  to 
the  unsteady  state  momentum  equation  (3.2)  are  plotted  in  Figs.  3.4 
and  3.5.  For  easier  comparison  with  previous  works  all  time  variations 
are  shown  plotted  as  a  function  of  t  by  using  the  inverse  of  the  trans¬ 
formation  (3.1). 

A  discontinuity  in  the  solution  was  also  encountered  in  the 
same  neighbourhood  of  t  as  in  the  results  of  Reference  [4]  and  [8]. 

An  attempt  has  been  made  in  these  works  to  discover  the  cause  of  the 
blow  up  condition  and  it  appears  that  the  discontinuity  occurs  along 
a  curve  defined  by  1  -  x  f£  =  0  (shown  in  Figs.  3.4  and  3.5)  which  is 
the  coefficient  of  2(Af ) '  in  equation  (2.35).  When  x  becomes  greater 
than  l/f£,  this  coefficient  changes  sign  from  positive  to  negative 
thereby  changing  the  partial  differential  equation  to  a  type  which, 
evidently,  is  no  longer  compatible  with  the  numerical  method  employed 
for  its  solution.  At  points  where  x  =  1/f'  in  equation  (2.30)  it  can 
be  readily  shown  that  the  mixed  derivative  9  f/9n9T  becomes  singular. 
The  very  fact  that  it  was  not  possible  to  obtain  a  continuous  solution 


■ 

- 

*  *  * 
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of  the  momentum  equation  to  its  final  steady  state  value  has  a  major 
influence  on  the  prime  objective  of  this  work,  namely,  to  obtain  a 
numerical  solution  to  the  unsteady  state  energy  equation  for  Prandtl 
numbers  different  from  unity. 

3.2  Numerical  Solution  of  the  Energy  Equation 

The  case  of  a  zero  pressure  gradient,  applicable  in  the  weak 
interaction  region,  allows  the  solution  of  the  momentum  equation  to 
be  obtained  independently  from  the  energy  equation.  For  the  special 
case  of  Pr  =  1  a  closed  form  solution  to  the  energy  equation  is  possible 
and  has  been  given  in  Reference  [4]  namely, 

g(n,T)  =  1  +^Me2[l  -  f'2(n,T)]  (3.21) 

for  zero  heat  transfer  at  the  wall  and 


g(n.x)  =  1  +  1^-He2f'(n,t)[l-f,(n,T)]+[gw-l][l-f,{n,-r)]  (3.22) 

for  non-zero  heat  transfer  at  the  wall.  Equation  (3.21)  satisfies  the 
differential  equation  (2.31)  and  also  the  boundary  conditions  (2.43) 
and  (2.45).  Similarly,  the  solution  (3.22)  satisfies  (2.31)  and  the 
boundary  conditions  (2.43)  and  (2.44).  It  is  implied  in  [4]  that  the 
initial  condition  corresponding  to  (2.46)  is  given  by 

=  1  +  ^  m22  [l-f ' (n,o)2] 


g(r|  ,0)  = 


(3.23) 
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for  zero  wall  heat  transfer,  or 

g(n,o)  =  =  l  +^^2f,fo,0)[l-f,(n,0)] 

+  [gw-i][l-f'(n,o)]  (3.24) 

for  the  case  of  heat  transfer  at  the  wall.  That  is  to  say  that  the 
enthalpy  jump  corresponding  to  the  M  velocity  jump  at  zero  time  can 
be  obtained  by  substituting  the  initial  condition  to  the  momentum 
equation  (2.42)  into  the  solutions  to  the  energy  equation  given  by 
(3.21)  and  (3.22). 

The  Prandtl  number  for  air  has  a  value  nearer  to  .72  than 
to  1  and  it  is  attempted  to  determine  what  effect  this  change  would 
have  on  the  solutions  (3.21)  and  (3.22).  The  steady  state  solution 
of  equation  (2.33)  for  arbitrary  Prandtl  numbers  is  discussed  in 
Section  3.3  and  it  is  well  known  that  reducing  Pr  from  1.0  to  .72  has 
a  significant  effect  on  the  enthalpy  distribution  in  the  boundary 
layer.  Figure  3.6  illustrates  this  for  a  free  stream  Mach  number  of 
10.  For  the  unsteady  energy  equation,  however,  a  closed  form  solution 
is  no  longer  possible  so  the  numerical  method  of  Section  3.1  has  been 
employed  to  solve  the  energy  equation  (3.3).  As  mentioned  previously, 
the  magnitudes  of  the  incremental  stream  function  Af  are  necessarily 
quite  small  and  also  small  in  relation  to  the  final  steady  state  dis¬ 
tribution  f2(n).  Similarly,  the  deviation  of  g(n,S)  from  the  steady 


f 


. 
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state  distribution  (n)  was  also  expected  to  be  small  in  magnitude. 
For  this  reason  the  author  has  further  rewritten  equations  (3.3)  in 
terms  of  Ag(n,£)  where 


Ag(ns£)  =  g2(n)  -  g(n,£)  (3.25) 

and  g2(n)  is  the  solution  to  equation  (2.33)  with  Mg  =  We  note 

that  the  motivation  behind  (3.25)  is  not  the  same  as  it  was  for  (2.34) 
since  (3.3)  is  already  a  linear  partial  differential  equation.  Thus 
in  terms  of  Ag  the  energy  equation  becomes 

{f-25(l-?)f}  |^-  2{l-5(l+f')}(l-5)  ff3- 
3n 

■  P?92'  ■  {f-250-5)hg^  -  (y-1)  Me2f"2  =  0  .  (3.26) 

The  finite  difference  approximations  used  are 

IIs-  =  27  {3Ag(m+l ,n)  -  4Ag(m,n)  +  Ag(m-1 ,n)  +  0(T2)  (3.27) 

|7=  2IT  {Ag(m+l,n+l)  -  Ag(m+l,n-l)}  +  0(H2)  (3.28) 


2 

\  (Ag(m+1 ,n+l)  -  2Ag(m+l ,n)  +  Ag(m+l,n-l)}  +  0(H2)  (3.29) 
3iT  H 


. 


and  the  computational  molecule  appears  as  sketched  in  Fig.  3.3 


Fig.  3.3  5-Point  Computational  Molecule 

Substituting  expressions  (3.27)  to  (3.29)  into  equation  (3.26)  and 
collecting  terms  one  obtains  with  Um  n  =  Ag(m,n) 


{PF+  T  -  {3H2[l-?(l+f)][l-5]  +  |I}  Um+1>n 

+  {  RT-  T  tMed-OfIJVun-l  =  -  4H2[l-S(l+f')][l-5]  Um>n 

2 

+  H2[l-S(l+f')][l-S]Um.ljn  +  ^  g“  +  H2T[f-25(l-S)f]g^ 

+  H2T(y-1)  He2f"2  .  (3. 


This  can  be  written  in  the  form  of  a  recursion 
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am+l  ,n  bm+l ,n+l  +  bm+l  ,n  bm+l  ,n  +  cm+l ,n  dm+l  ,n-l  dm+l  ,n  (3.31) 
with  Vl  ,n  =  T  {P?  +  H^r  '  5(l-Of]}  (3.32) 

bm+l,n  =  -  (3H2[l-5(l+f')][l-5]  -  fl  }  (3.33) 

cm+l,n  =  T{PF-H[f-50-5)f]}  (3.34) 

dm+l,n  =  H2i:i-5(1+f')][l-5][Um.ljn-4Umjn] 

g  •  ' 

+  H2T  {pf—  +  [f-25(l-5)f]g£  +  (Y-l)  Mg2f ' ,2}.  (3.35) 

In  an  analogous  manner  as  for  the  solution  to  the  momentum 
equation  each  (m+l)st  time  line  is  evaluated  according  to  the  recursion 
formula 


^rn+1  ,n  Rm+1  ,n  ^m+1  ,n+l  +  ^m+1  ,n 


(3.36) 


Now  from  (3.31) 


U 


m+1  ,n 


^m+1  ,n  “  cni+l  ,n  Um+l,n-l  ~  am+l  ,n  Um+1  ,n+l 

^m+1  ,n 


and  from  (3.36) 


Um+1 ,n-l  Rm+1 ,n-l  Um+1  ,n  +  Cm+1 ,n-l 
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so  that 


U 


_  ^m+1  ,n  Cm+1  ,nRm+1  ,n-lUm+1  ,n~Cm+l  ,n^m+l  ,n-l~am+l  ,n^m+l  ,n+l 
m+1 •"  "  bm+l  ,n 


or 


-a 


U 


=  [ 


m+1  ,n 


]  U, 


m+1  ,n  Lb  , ,  +c  ,1  R  ,  i  J  m+1  ,n+l 
*  m+1  ,n  m+1 ,n  m+1 ,n-l  * 


+ 


d  -r  C 
m+1  ,n  m+1  ,n  m+1  ,n-l 

^m+1 ,n+cm+l ,n^m+l ,n-l 


Upon  comparing  this  with  (3.36)  we  have 


m+1  ,n 

(3.37) 

m+1  ,n 

^m+1  ,n  +  cm+l  ,n 

^m+1  ,n-l 

d  —  c 

m+1  sn  m+1  ,n 

^m+1 ,n-l 

(3.38) 

m+1  ,n 

^m+1  .n  +  cm+l ,n 

Rm+1  ,n-l 

The  values  R  and  C  are  then  obtained  from  (3.37)  and  (3.38) 
m+1 ,n  m+l ,n 

by  marching  upwards  from  the  plate  towards  the  outer  edge  of  the 
boundary  layer  beginning  with  the  known  values  of  Rm+-j  i  and  C^+-j 
In  terms  of  Ag  boundary  condition  (2.44)  becomes 


Ag(O.e)  =  g2(0)  -  gw  =  0 


(3.39) 


and  from  (3.36) 


Um+1 ,1  =  Vl  ,1  Um+1 ,2  +  Cm+1 ,1 


(3.40) 


so  that  R  =  C  =  0  for  the  case  of  non-zero  heat  transfer  at 

m+1 ,1  m+1 ,1 


the  wall.  Boundary  condition  (2.45)  together  with  (3.25)  gives 
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Ag '  (0,0  =  g£(0)  -  g'(0,5)  =  0  (3.41) 

and  using  a  three-point  Newton  forward  difference 

Agm+l,l  "  2H  {_3Agm+l,l  +  4Agm+l ,2  ~  Agm+1 ,3} 

we  take  Agm+1 ^  =  Agm+^2  =  Agm+1  3  so  that  from  (3.40)  we  must  have 

*Vi+l  ,1  =  ^  anc*  ^m+l  1  =  g  ^or  case  zero  *iea^  transfer  at  the 
wall.  Having  obtained  all  n  and  Cm+-j  ^  we  apply  boundary  condition 
(2.43)  where  we  have 


Ag(ne»S)  =  0  . 


(3.42) 


The  values  of  Um+1  n  are  then  obtained  by  marching  down  towards  the 
plate  using  (3.36)  beginning  with  Um+-j  ^  =  0  . 


3.3  Determination  of  Stepsizes  H  and  T 

The  stepsizes  used  in  References  [4]  and  [8]  were  H  =  0.10 
and  T  =  0.05.  Rodkiewicz  and  Gupta  [8]  have  studied  the  effect  of  vary¬ 
ing  the  grid  size  in  the  solution  to  the  momentum  equation  and  found  that 
the  values  of  Af  changed  by  approximately  2.5%  when  T  and  H  were  chosen 
such  that  H  =  T  =  0.05.  A  similar  comparison  with  H  =  0.10  and  T  =  0.005 
resulted  in  a  5.8%  change  in  the  Af  values  i.e.  when  compared  with  the 
Af  values  obtained  with  the  original  grid  size  of  H  =  0.10  and  T  =  0.05. 


' 
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It  was  decided  in  Reference  [8]  to  keep  H  and  T  at  0.10  and  0.05  re¬ 
spectively,  in  order  to  preclude  large  round-off  errors  and  to  avoid 
large  computer  time  and  storage  requirements.  In  the  present  work, 
however,  it  was  found  that  the  total  error  produced  by  a  grid  size  with 
H  =  0.10  and  T  =  0.05  could  not  be  tolerated  and  furthermore  that  no 
large  round-off  errors  were  encountered  by  reducing  H  to  one-tenth  its 
original  value.  These  conclusions  are  based  on  the  following  calcu¬ 
lations  : 

Solutions  to  the  steady  state  energy  equation  (2.33)  are  well 
known  for  arbitrary  Prandtl  numbers.  Chapman  and  Rubesin  [20],  for 
example,  give  the  solution  first  advanced  by  Pohlhausen.  It  is 

g2(n)  =  1  +  r(n)  -  (y-1)  C,  s(n)  (3.43) 


where 


r(n)  =  2Pr 


oo 


[f£'(c)] 


Pr 


[f^'(e)]2'Pr  dedc 


(3.44) 


and 


s(n)  = 


oo 


c=n 


[f^'(?)]Pr  dc. 


(3.45) 


For  the  case  of  zero  wall  heat  transfer  we  take  C1  =  0  since  then 

dg2 


dr)  '  n=0 


=  0  so  that  for  the  insulated  plate  we  have  the  solution 


g2(n)  =  1  +  (*y~)  Me2  * 


(3.46) 


For  the  case  of  non-zero  wall  heat  transfer 


- 
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g2(0)  =  gw  =  1  +  ffp)  w\z  r(0)  -  (y-Dq  s(0) 


so  that 


ci  = 


1  +  t1?1)  Me2  r{0)  '  9 

(y-D  stO) 


w 


(3.47) 


Substituting  (3.47)  into  the  general  solution  (3.43)  we  obtain  for  the 
case  of  non-zero  wall  heat  transfer 

g2(n)  =  1  +  (-^)Me2  +  ^“^Me2  "  gw}  sjo}  (3-48) 

For  Pr  =  1  expression  (3.44)  reduces  to 


r(n)  =  1  -  [^(n)]' 


(3.49a) 


and 


s(n)  =  1  -  f£(n) 


(3.49b) 


so  that 


r(0)  =  s(0)  =  1  . 


Substitution  into  (3.46)  yields  for  zero  wall  heat  transfer  and  Pr  -  1 


g2(n)  =  1  +  ("^)Me2  ^  * 


(3.50) 


Similarly,  for  non-zero 


wall  heat  transfer  and  Pr  =  1  equation  (3.48) 
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becomes 


g2(n)  =  1  +  -  f£(n)]  +  [gw-l][l-f£(n)]  .  (3.51) 

We  note  immediately  that  for  Pr  =  1  the  solutions  (3.21)  and  (3.22)  to 
the  unsteady  energy  equations  are  exactly  of  the  same  form  as  the 
corresponding  solutions  (3.50)  and  (3.51)  to  the  steady  state  energy 
equation.  This  observation  has  been  further  investigated  and  leads 
to  a  plausible  alternative  solution  of  the  unsteady  energy  equation  for 
Prandtl  numbers  which  are  only  slightly  different  from  one.  Such  a 
solution  will  be  discussed  in  Section  3.5. 

Returning  to  the  problem  of  obtaining  an  optimum  stepsize  it 
was  decided  to  solve  the  steady  state  energy  equation  by  a  similar  re¬ 
currence  relation  as  used  in  the  solution  to  the  unsteady  momentum  and 
energy  equations.  Varying  the  stepsize  H  and  comparing  the  resulting 
numerical  solution  with  the  published  solution  (3.43)  would  be  an  indi¬ 
cation  of  the  magnitude  of  the  error  to  be  expected  in  the  numerical 
solution  of  equations  (3.2)  and  (3.26).  Figure  3.7  shows  the  error 
curves  obtained  in  the  numerical  solution  of  the  non-dimensional  surface 
equilibrium-enthalpy  for  a  free  stream  Mach  number  of  10,  y  =  1.4  and 
for  Prandtl  numbers  1  and  0.72.  For  Pr  =  0.72  the  "recovery  factor" 
r (0)  from  equation  (3.44)  is  equal  to  0.845  [20]  so  that  equation  (3.46) 
yields  g^ (0)  =  17.900.  Figure  3.7  shows  that  for  H  =  0.1  the  recurrence 
relation  method  resulted  in  an  error  of  approximately  4.0%  representing 


» 
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an  absolute  error  of  0.700.  For  Pr  =  1.0,  r(0)  =  1  so  that  g2(0)  = 
21.000  whereas  the  recurrence  relationship  gave  g2(0)  =  20.069  with 
H  =  0.1  which  is  an  error  of  about  4.5%.  Since  the  absolute  error 
in  these  cases  was  of  the  same  order  of  magnitude  as  Ag,  it  is  clear 
that  H  =  0.10  could  not  be  accepted.  Reducing  the  stepsize  H  to  one- 
tenth  of  this  value,  however,  resulted  in  an  error  of  less  than  0.5% 
in  both  cases.  Hence  it  was  decided  to  take  H  =  T  =  .01  for  the  nu¬ 
merical  solution  of  the  unsteady  momentum  and  energy  equations  (3.2) 
and  (3.26)  respectively.  Figure  3.7  shows,  furthermore,  that  the  ef¬ 
fects  of  round-off  errors  have  not  yet  significantly  come  into  play 
for  H  =  .01 . 

Finally,  it  should  be  pointed  out  that  the  actual  solutions 
to  the  steady  state  energy  equation  (3.43)  used  in  conjunction  with 
the  numerical  solution  of  equations  (3.26)  were  obtained  with  the 
application  of  Simpson's  rule  to  expressions  (3.44)  and  (3.45)  as 
supplied  by  a  scientific-subroutine-package  library  program  with  the 
IBM  360.  The  maximum  error  thus  introduced  was  at  most  0.20%  and  in 
the  majority  of  cases,  comparison  with  known  solutions  showed  a  dis¬ 
crepancy  of  less  than  0.1%. 

3.4  The  Assumption  g(n>0)  =  0  as  Initial  Condition  for  the  Unsteady 
Energy  Equation 

Examining  the  unsteady  energy  equation  (2.31)  it  is  seen 
that  we  require  two  boundary  conditions  in  space  and  one  in  time.  The 
space  boundary  conditions  (2.43)  and  (2.44)  or  (2.45)  apply  to  the  solu- 


. 
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ti on  of  (2.31)  for  any  arbitrary  Prandtl  number,  however,  this  is  not 
true  for  the  initial  condition  (2.46).  For  the  special  case  of  Pr  =  1, 
the  solution  given  by  (3.21)  or  (3.22)  was  obtained  in  Reference  [4] 
in  closed  form  without  any  specific  regard  to  the  initial  condition 
given  symbolically  by  (2.46).  In  fact,  the  initial  condition  is  implied 
from  the  solution  and  is  given  by  either  (3.23)  or  (3.24).  For  Pr 
different  from  unity,  a  closed  form  solution  to  the  unsteady  energy 
equation  cannot  be  found  and  hence  it  becomes  necessary  to  resort  to 
numerical  techniques.  The  numerical  method  has  already  been  described 
in  Section  3.2  and  we  recall  that  the  step  by  step  procedure,  as  with 
any  numerical  technique,  requires  a  given  initial  time  line  before  it 
is  able  to  commence. 

It  is  evident  that  the  solution  g(n,T)  obtained  from  (2.31) 
must  eventually  coincide  with  some  final  steady  state  solution  g2(n) 
in  the  limit  as  x  ■*  °°  (?  -*  1)  and  this  would  satisfy  the  requirement 
of  one  boundary  condition  in  time.  Intuition  indicates  that  the  ap¬ 
proach  to  the  final  steady  state  would  be  asymptotical  so  that  a  method 
similar  to  Reference  [25 J  might  be  applied  where  an  original  guess  at 
the  initial  g(n50)  distribution  is  successively  modified  until  the 
boundary  condition  at  infinity  is  satisfied  to  within  a  predetermined 
tolerance  and  sufficiently  large  enough  value  of  x.  Unfortunately 
any  thought  in  this  direction  is  already  ruled  out  since,  as  was 
pointed  out  before,  the  solution  to  the  unsteady  momentum  equation 
is  not  available  for  large  values  of  x.  Thus,  lacking  any  experimental 
evidence,  it  becomes  necessary  to  fall  back  on  pnysical  and  mathematical 
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reasoning  to  arrive  at  a  reasonable  approximation  to  the  initial  en¬ 
thalpy  distribution  g(n,0)  corresponding  to  a  sudden  increase  in 
velocity  in  the  boundary  layer. 

Now  any  solution  g(n,x)  to  equation  (2.31)  must  necessarily 
satisfy  this  equation  for  all  times  and  in  particular 

2Pr  |3- -  {g"  +  Prfg'  +  Pr(y-l)  Me2f2}  =  0  (3.52) 

which  is  the  unsteady  energy  equation  (2.31)  specialized  to  i  =  0.  Be¬ 
fore  any  jump  in  velocity  takes  place,  the  enthalpy  distribution  in 
the  boundary  layer  is  given  by  the  solution  of  the  initial  steady 
state,  namely,  the  solution  of 

9q'  +  pr  <o  9q  +  My-D  Me02  fy2  =  0  (3.53) 

where  subscript  "0"  means  symbolically  t  <  0.  Intuitively  one  might 

think  that  the  enthalpy  jump  corresponding  to  a  sudden  increase  in 

velocity  would  simply  be  given  by  the  solution  of  (3.53)  with  the 

initial  values  of  the  momentum  equation  substituted  in  place  of  the 

Blasius  values.  If  this  were  indeed  so,  it  would  imply  that  3g/9x(ri ,0) 

g(n ,0)  =  0  in  expression  (3.52).  In  view  of  definition  (3.25)  for  Ag 

we  have  9Ag/9x  =  9Ag/9£  =  -  9g/9x  at  x  =  0.  When  the  solutions  (3.21) 

and  (3.22)  for  Pr  =  1  are  plotted  in  terms  of  Ag  as  a  function  of  x 

as  in  Fig.  3.8  for  the  insulated  wall  and  also  for  non-zero  wall  heat 

transfer  (case  q  =  0),  it  is  observed  that  in  both  cases  9Ag/9x(n,0) 
v  3w 


■ 


■ 
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is  virtually  zero  even  for  very  small  values  of  q.  We  recall  that 
one  of  the  basic  hypersonic  assumptions  (2.24)  already  made  states 
that  3h  /8t  =  0  so  that  in  the  free  stream  3g/3x  =  0  is  satisfied 

c 

exactly.  We  have  also  seen  that  for  the  case  of  Pr  =  1  the  g-distribu- 
tion  is  obtained  at  t  =  0  simply  by  specializing  the  solution  to  the 
energy  equation  to  zero  time,  that  is,  the  initial  values  of  the  mo¬ 
mentum  equation  are  substituted  in  the  solution  to  obtain  the  initial 
conditions  (3.23)  and  (3.24).  Figure  3.9  shows  the  results  obtained 
when  the  initial  condition  (3.24)  with  gw  =21  is  used  in  the  numerical 
solution  of  the  energy  equation  with  Pr  =  1  up  to  x  =  1 .  This  is 
compared  with  the  closed  form  solution  (3.22),  in  terms  of  Ag,  and  it 
is  clear  that  the  two  solutions  agree  almost  exactly. 

For  an  arbitrary  Prandtl  number  there  is,  of  course,  no 
closed  form  solution  and  an  initial  g(n,0)  distribution  must  be  found 
in  a  different  way.  Rewriting  the  exact  equation  (3.52)  at  zero  time 
in  the  form 


g"  +  Prfg’  =  Pr[2g  -  (y-1 )  (M0f '  * )2]  (3.54) 

and  assuming  for  the  moment  that  the  correct  g(q,0)  was  known  from 
some  other  source,  the  entire  right-hand  side  of  (3.54)  could  then  be 
yj gwed  as  the  inhomogeneous  part  of  a  second  order  ordinary  differential 
equation  for  g.  Now  for  Pr  =  1  it  is  formally  possible  to  obtain  the 
range  of  magnitude  for  g(o,0);  we  have  already  shown  an  example,  by 
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plotting  the  solutions  in  Fig.  3.8,  where  g(n,0)  =  -  Ag(n,0)  is  small 

for  all  n-  An  exact  calculation  for  a  particular  case  of  q  =21, 

Pr  =  1,  y  =  and  Mg  =  10,  shows  that  g(o,0)  reaches  its  maximum 

value  of  only  0.9  very  close  to  the  wall  (n  =  .04),  diminishes  rapidly 

to  zero  between  q  =  .43  and  .44  and  attains  a  magnitude  no  greater 

than  .056  as  it  approaches  zero  at  n  =  n  .  Thus,  for  Pr  =  1 ,  the 

0 

largest  magnitude  of  2g(n,0)  is  1.8.  On  the  other  hand,  f"(n,0)  = 
f 2 ' ( n )  +  Af,,(n,0)  is  also  larger  in  magnitude  near  the  plate  and  also 
tending  to  zero  at  n  =  ri  •  In  particular,  at  n  =  .04  f 1 '  =  .4649  so 

G 

o 

that  the  term  (y-l)(M  f11)  is  of  magnitude  8.64  or  approximately 

five  times  as  great.  In  an  order  of  magnitude  analysis  this  would 

•  p 

not  quite  justify  neglect  of  the  term  2g  relative  to  (y-l)(M  f 1 ' )  . 

e 

However,  only  slightly  further  up  from  the  wall  and  out  to  the 
boundary  layer  edge  the  two  terms  compare  more  like  45  to  1  which 
clearly  justifies  neglect  of  2g(n,0)  in  equation  (3.54)  for  the  major 
portion  of  the  boundary  layer  excluding  only  the  region  very  close 
to  the  plate  surface. 

Considering  the  observations  made  above  for  the  case  of 
Pr  =  1,  it  is  speculated  that  for  arbitrary  Prandtl  numbers  different 
from  unity,  the  required  g-distribution  at  zero  time  might  be  obtained 
from  the  solution  of 

gj'(n,0)  +  Pr  f, (n,0)  gj(n.O)  +  Pr(y-l)  fj,2(n,0)  =  0  (3.55) 


where  subscript  "1"  means  symbolically  t  =  0.  Thus  we  are  substituting 
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the  initial  values  of  the  momentum  equation  into  equation  (3.52)  and 
make  the  assumption  that  g(n,0)  =  0  for  arbitrary  values  of  Pr. 

From  physical  considerations,  the  assumption  that  g(n,0)  =  0 
is  further  reinforced  when  we  recall  Section  2.5  where,  for  the 
momentum  equation,  it  was  argued  that  at  the  instant  of  change  of 
the  free  stream  velocity  the  governing  equation  in  the  major  portion 
of  the  boundary  layer  was 


3u  =  _  1 
3t  ’  p  3x 


(3.56) 


Since  we  are  dealing  with  the  region  of  zero  pressure  gradient  we 
have,  for  the  momentum  equation,  3u/9t  =  0  at  zero  time.  If  we  apply 
the  same  reasoning  to  the  energy  equation  (2.4)  we  would  neglect  all 
v.iscous  and  convective  terms  and  again  for  constant  pressure  the 
governing  energy  equation  at  zero  time  would  then  be 


fib0  (3-57) 

which  is  equivalent  to  saying  3g/3T(n>0)  =  0  in  equation  (2.31). 

3.5  Alternative  Solution  of  the  Unsteady  Energy  Equation 

It  was  observed  in  Section  3.3  that  the  form  of  the  solutions 
for  the  steady  state  enthalpy  distributions  (3.50)  and  (3.51)  are 
identical  with  the  corresponding  unsteady  distributions  (3.21)  and 
(3.22)  when  the  Prandtl  number  equals  unity.  That  is,  the  time  vari- 
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able  t  in  the  solutions  appears  only  as  an  independent  variable  in  the 
function  f'(n>T).  When  the  Prandtl  number  does  not  differ  greatly 
from  unity  as  in  the  case  of  gases,  N.  Curie  [26]  published  an  empirical 
steady  state  enthalpy  distribution  (due  partly  to  R.J.  Monaghan)  which 
satisfies  the  boundary  conditions  both  at  the  plate  surface  and  at 
the  edge  of  the  boundary  layer.  Thus  for  0.725  £  Pr  <_  1.25  and  with 
the  notation  used  here  the  empirical  distributions  become 

g  =  1  +  [Pr1/2  ^M^HMf')3]  -  £LMe2Pr[(f')2-(f1)3]  (3.58) 

for  the  adiabatic  wall,  and 

g  =  l+[gw-l][l-(f')3]  +  Pr1/3[1+Pr1/2  ^  Me2-gw][f'-(f')3] 

(3.59) 

-  ^Me2  Pr[(f')2  -  (f')3] 
for  the  isothermal  wall. 

Figure  3.10  shows  a  comparison  of  the  enthalpy  distributions 
obtained  from  the  empirical  formulae  and  from  the  exact  solution  given 
by  equation  (3.43)  for  the  case  of  Pr  =  0.72  and  M^  10.  Agreement 
is  seen  to  be  fairly  good  throughout  the  entire  boundary  layer  and  al¬ 
most  exact  agreement  is  obtained  in  the  more  important  region  close  to 

the  plate  surface. 

To  justify  the  use  of  expressions  (3.58)  and  (3.59)  as 
approximate  solutions  of  the  unsteady  energy  equation  (3.3)  it  is 
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necessary  to  show  that  the  degree  of  satisfaction  of  (3.3),  even  though 
not  exact,  is  at  least  constant  throughout  the  time  domain  of  integra¬ 
tion.  That  equation  (3.3)  would  not  be  satisfied  exactly  by  such  an 
empirical  distribution  is  evident  from  Fig.  3.10  although  better  satis¬ 
faction  would  be  expected  nearer  the  plate.  If  we  take  for  example  the 
case  of  zero  wall  heat  transfer  with  Pr  =  0.72,  =  10,  and  x  =  0.25, 

p  p 

evaluate  from  (3.58)  the  partial  derivatives  3g/3q,  3  g/3n  »  9g/9£, 
substitute  these  into  equation  (3.3),  denote  by  D  the  degree  of  satis¬ 
faction  (i.e.  the  deviation  from  zero)  of  (3.3)  and  compare  this 
quantity  with  the  highest  order  partial  derivative,  we  obtain  the  re¬ 
sults  shown  in  Fig.  3.11.  It  is  seen  that  satisfaction  is  not  exact 
(D  =  0)  as  expected  but  the  degree  of  satisfaction  is  about  the  same 
for  the  unsteady  energy  equation  as  for  the  initial  steady  state. 

Figure  3.12  shows  the  quantity  D  plotted  against  x  and  it  is  observed 
to  be  practically  constant  in  the  range  of  x  from  0  to  1  for  each  n 
with  the  degree  of  satisfaction  improving  as  we  approach  the  wall.  A 
similar  result  was  observed  for  the  isothermal  wall  using  expression 
(3.59)  for  g(n,£). 

In  concluding  this  Chapter,  it  is  worth  pointing  out  that 
aside  from  verifying  that  (3.24)  is  the  correct  g(n,0)  distribution 
for  the  solution  (3.22)  Fig.  3.9  also  serves  to  show  that  the  error  in 
the  numerical  procedure  used  to  solve  the  unsteady  equations  is  quite 
small,  as  was  expected  from  the  results  of  Section  3.3. 
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CHAPTER  IV 

TRANSIENT  HEAT  TRANSFER,  SKIN  FRICTION 
AND  INDUCED  PRESSURE 

4.1  Unsteady  Heat  Transfer  Coefficient 

The  heat  transfer  at  the  plate  is  given  by 


k 


y=0 


or  using  (2.7)  and  (2.17) 
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(4.1) 


(4.2) 


In  high  speed  flow  it  has  been  found  experimentally  that 
the  direction  of  heat  flow  at  the  wall  does  not  depend  on  the  dif¬ 
ference  between  the  wall  temperature  and  the  free  stream  temperature 
as  in  sub-sonic  flow  but  rather  on  the  difference  between  the  wall 
temperature  and  the  adiabatic  wall  temperature.  Thus  the  surface  con¬ 
vective  conductance  is  defined  by 


a 


raw' 


(4.3) 


where  a  is  time  dependent  for  the  unsteady  problem,  or  in  other  words. 
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a  represents  the  instantaneous  local  heat-transfer  coefficient  [27]. 
T*raw  is  the  reference  adiabatic  wall  temperature  which,  for  consistency 
with  the  steady  state,  is  taken  as  the  final  steady  state  adiabatic 
wall  temperature.  Using  (4.2)  and  the  definitions  of  the  Prandtl 
number  and  local  Reynolds  number 


(4.4a) 


and 


V 

v 


(4.4b) 


we  obtain 


(/§)  9'(0,t)  =  (/  Rex)  Pr  St(graw-gw) 
where  the  Stanton  number  has  also  been  employed  with  the  definition 
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In  terms  of  the  Nusselt  number  Nu  =  St  Re  Pr  the  local,  non-dimensional, 
unsteady  heat  transfer  coefficient  may  then  be  written  as 


Nu. 
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(4.6) 
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4.2  Unsteady  Skin  Friction  Coefficient 

The  skin  friction  coefficient  Cf  is  defined  as 


C.  = 


w 


Z-  pref  Ue 


(4.7) 


where  the  wall  shear  stress  t  is  qiven  by 

w  3  J 


y=0 


(4.8) 


Again  using  (2.7)  and  (2.17)  we  obtain 
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(4.9) 


so  that  (4.7)  becomes 


Re. 


Cf  /  =  f 1 1  (0 ,x)  . 


(4.10) 


4.3  The  Vleak-Interaction  Unsteady  Induced  Pressure 

For  the  acoustic  approximation  to  the  induced  pressure  we 
have  from  [12] 


pe  aeve 


(4.11) 


From  [15]  the  expression  for  with  Ug  a  constant  is 
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v  =  PA  =  3A  +  8A 
e  Dt  at  e  8x 


Also  from  [15]  the  expression  for  the  displacement  thickness 
unsteady  boundary  layer  is 
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For  our  case  U  and  p  are  constants  so  that  (4.13)  becomes 
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where 
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and 
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In  terms  of  the  transformed  variables  used  these  expressions 

in  turn 
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of  the 
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become 
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6  = 
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(4.18) 


Combining  (4.11),  (4.12)  and  (4.14)  we  obtain 


=  {!^+u  36* 

PQ  P„  at  +  ue  3x 
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or 
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Differentiating  (4.17)  with  respect  to  x  and  (4.18)  with  respect 

to  t 


§f  =(/  ^■){-i-+{/x)ff 

3X  Ue  2/x  3X 


(4.21) 
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where 
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I  =  I (t)  =  f  (g-f')dn 


(4.23) 
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Introducing  the  fundamental  hypersonic  interaction  parameter  x  [1,3] 
defined  as 
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then 
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and 
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Substituting  (4.26)  and  (4.27)  into  expression  (4.20)  we  obtain 
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which  for  the  final  steady  state  reduces  to 


(4.28) 
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CHAPTER  V 

DISCUSSION  AND  CONCLUSIONS 
5.1  Discussion  of  Results 

Keeping  in  mind  the  small  magnitudes  of  our  Ag  function  de¬ 
fined  by  (3.25),  it  had  to  be  ensured  that  the  errors  involved  in  the 
numerical  solutions  to  both  g2(n)  and  g-j(n»0)  as  well  as  those  in  the 
solution  of  the  unsteady  energy  equation  (3.26)  are  by  far  below  the 
scale  of  Ag.  We  have  already  seen  from  Fig.  3.9  that,  using  the  closed 
form  initial  condition,  the  error  in  the  solution  of  (3.26)  for  Pr  =  1 
is  virtually  negligible.  In  Section  3.3  it  was  shown  that  the  errors 
in  g2(u)  were,  relatively  speaking,  insignificant  as  well.  Although 
equation  (3.55)  is  of  the  same  form  as  equation  (2.33)  the  solution 
given  by  (3.43)  does  not  apply  since  f-j(n,0)  does  not  satisfy 
f .j '  ■  +  f-jf-j'  =  0.  However,  the  recurrence  relationship  used  to  deter¬ 
mine  the  stepsize  H  in  Section  3.3  was  used  as  a  means  of  calculating 
the  g-j(ns0)  distribution  from  equation  (3.55)  with  a  stepsize  H  =  0.01. 
In  addition  the  g-j(n>0)  distribution  associated  with  (3.55)  was  also 
obtained  from  a  general  subroutine  called  DLBVP ,  which  is  supplied 
with  the  Scientific  Subroutine  Package  for  the  IBM  360  [28],  and 
solves  a  general  linear  boundary  value  problem.  The  results  obtained 
with  the  recurrence  method  checked  almost  perfectly  with  those  using 
this  subroutine.  The  added  advantage  of  subroutine  DLBVP  is  that  it 
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permits  the  subdivision  of  the  initial  stepsize  H  =  0.01  as  much  as 
ten  times  in  order  to  attain  an  almost  arbitrary  accuracy  specified 
by  the  user.  Thus  we  can  be  reasonably  confident  that  any  discre¬ 
pancies  observed  in  subsequent  results  are  most  likely  to  be  due  to 
our  assumptions  or  approximations  and  not  due  to  the  numerical  methods. 

Recalling  that  the  underlying  concern  of  this  work  was  to 
obtain  an  enthalpy  "jump"  in  the  boundary  layer  corresponding  to  a 
certain  jump  in  the  free  stream  velocity,  the  assumption  g(n,0)  =  0 
should  really  be  scrutinized  in  the  following  light:  One  of  the 
fundamental  characteristics  of  hypersonic  flow  is  that  relatively 
small  perturbations  in  the  velocity  are  accompanied  by  relatively 
large  perturbations  in  the  density,  pressure,  temperature,  and  speed 
of  sound  [29].  Hence  in  hypersonic  flow  the  Mach  number  changed  pri¬ 
marily  as  a  result  of  the  change  in  the  speed  of  sound.  Let  us  con¬ 
sider  as  an  example  a  flow  for  which  the  final  free  stream  Mach  number 
Mg2  has  a  value  of  10.  Then  it  may  be  shown  that  before  the  U  jump 
in  velocity  occurred,  the  initial  free  stream  Mach  number  Me0  was 
approximately  8.45  for  a  gas  with  y  =  1.4.  For  the  case  of  Pr  =  1 

and  g  =  0  say,  we  can  calculate  from  equation  (3.53)  the  initial 
w 

steady  state  enthalpy  distribution  gQ (n )  corresponding  to  a  Mach 
number  of  8.45  by  using  solution  (3.22).  From  the  closed  form  solu¬ 
tion,  equation  (3.24),  we  can  also  calculate  the  g-j  (rj ,0)  distribu¬ 
tion  corresponding  to  Me-j  =  M^  =  10.  Subtracting  gQ (n )  from  g-|(n»0) 
will  yield  the  corresponding  jump  in  enthalpy  for  the  case  of  Pr  =  1. 
The  solid  curve  in  Fig.  5.1  shows  the  resulting  enthalpy  jump.  If 
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we  now  apply  our  assumption  of  g(n,0)  =  0  and  hence  calculate  the 
g1 (n .0)  distribution  from  equation  (3.55)  we  obtain  an  enthalpy  jump 
shov/n  plotted  by  the  broken  curve  in  Fig.  5.1.  The  corresponding 
enthalpy  jumps  for  the  adiabatic  wall  is  shown  in  Fig.  5.2.  The 
first  thing  to  notice  is  that  the  enthalpy  jump  is  certainly  not  a 
constant  as  might  perhaps  be  speculated  from  the  discussion  in  Section 
3.4.  The  second  observation  is,  of  course,  that  our  assumption  re¬ 
sults  in  an  enthalpy  jump  which  closely  resembles  that  predicted  by 
the  closed  form  solution,  at  least  within  errors  that  could  be  tolerated. 
From  Fig.  5.1  we  determine  a  maximum  discrepancy  of  about  7.5%  and 
from  Fig.  5.2  about  9.5%.  It  is  most  likely  that  any  experimental 
determination  of  the  unsteady  enthalpy  profiles  would  contain  errors 
of  magnitude  at  least  of  the  same  order  as  depicted  in  Figs.  5.1  and 
5.2. 

In  order  to  discuss  the  effect  of  the  deviation  of  the 
Prandtl  number  from  unity,  solutions  for  Pr  =  1  should  be  obtained 
by  the  application  of  the  same  assumptions  or  approximations  as  for  the 
solutions  of  Pr  different  from  one.  That  is,  in  Section  3.4  all 
comparison  should  be  made  on  the  basis  that  g (n , 0 )  =  0  and  in  Section 
3.5  comparison  should  be  made  with  the  closed  form  solutions  of 
Reference  [4]  since  (3.58)  and  (3.59)  reduce  to  (3.21)  and  (3.22) 
respectively  when  Pr  =  1 . 

Figures  5.3  and  5.4  illustrate  the  effect  of  reducing  the 
Prandtl  number  from  1  to  0.72  upon  the  assumption  that  g(n50)  -  0. 
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The  solutions  of  (3.26)  here  for  the  adiabatic  wall  and  the  cold  v/all 
(9W  _  0)  have  been  plotted  in  the  form  Ag (ri ,t )  =  g(n>T)  -  g-j(ns0)  since 
it  is  not  entirely  obvious  that  the  solution  g2(n)  of  (2.33)  is  the 
corresponding  final  steady  state  when  the  initial  condition  g1(n,0) 
is  given  by  the  solution  of  (3.55).  The  initial  steady  state  gQ(n) 
is  a  natural  reference  state  for  this  solution  but  presentation  of  the 
results  in  the  form  Ag(n,x)  =  g(n,x)  -  gQ(n)  did  not  lend  itself  very 
well  in  displaying  the  transient  behaviour  of  the  solution.  The  in¬ 
stantaneous  local  Nusselt  number  given  by  expression  (4.6)  is  repre¬ 
sented  in  Fig.  5.5  with  M  =  10  and  for  Pr  =  1.0  and  0.72.  That  the 

-  c 

enthalpy  gradient  at  the  wall  g^O)  is  less  for  Pr  =  .72  than  for 
Pr  =  1.0  is  evident  from  inspection  of  Fig.  3.6,  however,  the  reference 
enthalpies  g^aw  differ  also  for  Pr  =  .72  (gfaw  =  17.9,  =  10)  and 

for  Pr  =  1.0  (g  all  =  21,  M  =  10).  Thus  it  might  have  been  antici- 
pated  that  (4.6)  is  independent  of  Prandtl  number.  A  net  reduction 
of  the  Nusselt  number,  however,  makes  physical  sense  since  the  dynamic 
viscosity  y  is  smaller  for  a  gas  with  a  lower  Prandtl  number  which  in 
turn  results  in  a  lower  frictional  heating  of  the  gas  and  hence  a 
lower  rate  of  heat  transfer  to  the  wall  when  the  wall  enthalpy  gw  is 
less  than  the  adiabatic  wall  enthalpy.  We  note  also  in  Fig.  5.5  that 
the  curves  differ  slightly  with  different  values  of  gw  and  that  the 
coefficient  rises  rapidly  to  a  more  or  less  constant  value.  In  view 
of  Fig.  3.4  the  point  t  =  n  =  0  must  be  excluded  for  which  reason  the 
curves  on  Fig.  5.5  have  not  been  extended  to  i  =  0. 

The  transient  weak-interaction  induced  pressure  given  by 
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expression  (4.28)  is  shown  plotted  in  Fig.  5.6  for  Pr  =  1  and  .72  with 
X  -  4.5  and  -  10.  The  value  of  x  =  4.5  has  been  chosen  merely  for 
comparison  of  results  given  in  Reference  [4];  for  the  weak  interaction 
region  the  interaction  parameter  x  assumes  a  magnitude  not  much  greater 
than  one  [1,3].  Again  it  appears  as  if  the  transient  portions  of  the 
curves  may,  for  all  practical  purposes,  be  neglected,  at  least  in 
the  event  of  only  a  small  velocity  impulse  of  1%.  Strickly  speaking 
the  points  x  =  0  should  again  be  disregarded  and  it  has,  furthermore, 
been  pointed  out  [12]  that  the  acoustic  approximation  (4.11)  becomes 
invalid  at  x  =  0. 

The  results  from  the  empirical  distributions  of  Section  3.5 

are  shown  on  Figures  5.7  and  5.8.  They  are  plotted  in  the  form 

Ag(n,x)  =  g^ri)  -  g(o}x)  this  time  since  the  solutions  are  made  to 

satisfy  the  necessary  time  boundary  condition  at  infinity.  The  solid 

curves  for  Pr  =  1  have  also  been  obtained  from  (3.58)  and  (3.59)  but 

they  are  in  fact  coincident  with  the  closed  form  solutions  given  in 

Reference  [4].  An  interesting  outcome  of  the  numerical  solutions  for 

g(n,x)  pertains  to  the  case  of  an  insulated  wall.  The  boundary  condition 

to  be  satisfied  at  the  wall  is  3g/3n (0 ,x ) |  =  0  i.e.  (2.45),  but  also, 

n=0 

because  there  is  no  heat  transfer  to  the  wall  and  no  slip  velocity, 
the  wall  enthalpy  must  remain  constant.  The  velocity  of  the  layer  of 
fluid  immediately  adjacent  to  the  wall,  however,  is  varying  as  the 
boundary  layer  adjusts  itself  to  the  new  steady  state  so  that  the 
temperature  in  this  layer  is  also  expected  to  vary  during  this  period 
of  readjustment.  The  result  is  a  temperature  jump  at  the  wall  which 
is  brought  out  by  the  computer  program  because  the  finite  difference 
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representati on  of  the  derivative  at  the  wall  relates  the  wall  temper¬ 
ature  to  the  temperature  in  the  first  layer  of  fluid.  Figure  5.3  de¬ 
picts  such  a  jump  whereas  Fig.  5.7  does  not.  We  also  note  that  the 
associated  boundary  condition,  namely  3Ag/3n(0,-r)  =  0  is  satisfied 
for  all  time  curves  in  our  solution  assuming  g(n,0)  =  0  whereas  this 
boundary  condition  is  not  satisfied  in  Fig.  5.7  for  t  =  0.  Thus,  some 
doubt  is  cast  on  the  validity  of  the  closed  form  solution  and  it  must 
be  emphasized  that  without  experimental  evidence  it  is  difficult  to 
pass  judgement  as  to  which  of  the  two  solutions  predicts  more  closely 
the  actual  state  of  affairs. 

The  heat  transfer  expression  (4.6)  from  the  empirical  solu¬ 
tion  is  shown  in  Fig.  5.9.  Here  also  a  rapid  advance  to  a  constant  value 
is  observed  although  high  values  are  predicted  for  small  t  in  contrast 
with  Fig.  5.5.  Also  the  two  curves  in  Fig.  5.9  represent  all  values 

of  M  and  g  which  was  not  the  case  of  Fig.  5.5.  That  the  heat 
e  w 

transfer  data  is  obtainable  from  a  single  curve  for  Pr  =  1  is  the  re¬ 
sult  of  Reynolds  analogy  and  it  can  be  shown  that  the  solid  curve  in 
Fig.  5.9  may  be  obtained  from  the  shear  stress  parameter  in  Fig.  5.13. 

For  Pr  =  0.72,  Reynolds  analogy  is  no  longer  valid  but  the  fact  that 
we  have  based  the  Mussel t  number  on  the  difference  T^  -  T^aw  renders 
the  Mussel t  number  as  some  function  f(Pr,Rex)  i.e.  explicit  dependence 
only  on  the  Reynolds  and  Prandtl  numbers.  Conversely,  basing  the 

Mussel t  number  on  the  difference  T  -  T  would  lead  to  some  function 

w  e  (  -1)  2 

f(Pr,Re  ,E)  i.e.  explicit  dependence  on  the  Eckert  number  E  ='(g  -{)  He  [2]. 
x  ■  w 
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The  weak- interaction  induced  pressure  is  shown  in  Fig.  5.10 

again  with  M0  =  10  and  x  =  4.5.  These  curves  are  very  similar  to  those 

in  Fig.  5.6  except  that  for  smaller  g  the  transition  to  a  constant 

w 

value  is  sligntly  more  gradual .  One  major  difference,  however,  is 
that  for  gw  ^_10  the  induced  pressure  with  the  assumption  g(n,0)  =  0 
for  Pr  =  .72  is  larger  than  for  Pr  =  1.0  whereas  in  Fig.  5.10,  using 
the  empirical  distribution,  the  induced  pressure  is  smaller  for  Pr  = 

.72  for  all  values  of  gw  up  to  its  equivalent  adiabatic  value.  If  we 
calculate  a  steady  state  induced  pressure  from  equation  (4.29)  using 
the  solution  g2(n)  of  equation  (2.33)  as  the  final  steady  state  en¬ 
thalpy  distribution  and  plot  the  results  as  a  function  of  gw  for  both 
Pr  =  1  and  Pr  =  0.72  the  curves  shown  in  Fig.  5.11  are  obtained.  We 
observe  that  the  induced  pressure  for  the  steady  state  case  and  for 
Pr  =  .72  does  indeed  exceed  that  for  Pr  =  1  for  values  of  gw  greater 
than  approximately  10.  Vie  note  also  that  the  solution  of  (2.33)  given 
by  (3.43)  is  exact  and  is  free  from  the  kind  of  assumptions  and  approxi¬ 
mations  made  for  the  unsteady  case.  Thus  Fig.  5.6  appears  to  give  to 
more  correct  results  although  it  is  possible  in  Fig.  5.10  that  the 
Pr  =  .72  curves  for  gw  >  10  might  eventually  overtake  the  Pr  =  1.0 
curves  since  they  have  not  yet  exactly  leveled  off  at  x  =  1.  Experi¬ 
mental  evidence  would  again  be  helpful  here. 

Figure  5.12  shows  the  transient  induced  pressure  as  a  function 
of  the  interaction  parameter  x  for  gw  =  0  and  for  different  free  stream 
Mach  numbers  as  obtained  from  the  empirical  distribution.  The  Mach 
number  is  seen  to  have  only  a  very  small  effect  when  increased  from 
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5  to  10,  which  is  a  rather  surprising  result.  Finally  Fig.  5.13  shows 
the  wall  shear  stress  parameter  of  Section  4.2  which  differs  only 
slightly  from  the  curve  given  in  Reference  [4]  since  it  was  calculated 
here  using  the  10-point  computational  molecule. 

5.2  Conclusions 

In  the  steady  state  case,  a  reduction  of  Prandtl  number  from 
1  to  0.72  has  the  effect,  according  to  Fig.  3.6,  of  reducing  the 
maximum  dimensionless  enthalpy  by  approximately  19.5%  for  the  cold 
wall  (g  =  0)  and  approximately  15%  for  the  adiabatic  wall.  According 
to  the  assumption  g(n>0)  =  0,  the  corresponding  reduction  in  the 
transient  Ag  distributions  is  on  the  average  approximately  8%  for  the 
adiabatic  wall  and  15%  for  the  cold  wall.  With  the  empirical  dis¬ 
tributions  we  observe  an  approximate  average  reduction  of  16%  for  the 
adiabatic  wall  (Fig.  5.7)  and  22%  for  the  cold  wall  (Fig.  5.8). 

When  viewed  in  this  manner  the  Prandtl  number  effect  is  about  the 
same  for  the  unsteady  solutions  as  for  the  steady  state  case  except 
perhaps  for  the  results  of  Fig.  5.3. 

The  effect  on  the  local  heat  transfer  expression  (4.6)  in 
accordance  with  the  assumption  g(n,0)  =  0  (Fig.  5.5)  is  a  reduction 
by  about  11%  for  gw  =  0,  11.3%  for  gw  =  10  and  12.2%  for  gw  =  15  at 
t  =  1.  Similarly  from  Fig.  5.9  using  the  empirical  distribution, 
the  reduction  due  to  a  change  in  Pr  from  1  to  .72  is  10.3%  at  t  1 
which  applies  for  all  M  and  g  below  the  equivalent  adiabatic  value. 
The  assumption  g(n,0)  =  0  predicts  a  fairly  sudden  increase  in  heat 
transfer  from  some  lower  value  at  t  =  0  to  a  larger,  more  or  less 
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constant  value.  The  empirical  distribution  results  of  Fig.  5.9  for 
Pr  =  .72  on  the  other  hand  follows  the  form  of  the  closed  form  solu¬ 
tion  for  Pr  =  1  which  predicts  a  slightly  more  gradual  reduction  in 
heat  transfer  from  some  larger  value  at  i  =  0  to  a  more  or  less  con-' 
stant  value  at  later  times. 

The  induced  pressures  from  Figs.  5.6  and  5.10  behave  quite 
similar  and  it  appears  that  the  transient  behaviour  may  be  neglected 
since  a  fairly  constant  value  is  reached  almost  instantaneously. 

The  results  with  g(q,0)  =  0  predict  a  behaviour  with  g  which  is  more 
consistent  with  that  observed  in  the  final  steady  state,  Fig.  5.11, 
namely,  that  a  larger  induced  pressure  is  observed  with  a  smaller 
Prandtl  number  for  wall  enthalpy  ratios  g  greater  than  about  10  but 
less  than  their  equivalent  adiabatic  wall  values.  With  the  assumption 
of  9 (n ,0 )  =  0  a  reduction  in  the  induced  pressure  according  to  Fig. 

5.6  is  observed  as  follows:  8.4%  for  the  adiabatic  wall  and  about  15% 
for  the  cold  wall.  From  Fig.  5.10  we  observe  a  reduction  of  12.9% 
for  the  adiabatic  wall  and  22.6%  for  the  cold  wall. 

Finally  we  conclude  that  the  Mach  number  effect  on  the 
induced  pressure  as  seen  from  Fig.  5.12  is  virtually  negligible  com¬ 
pared  with  the  effect  of  a  small  change  in  Prandtl  number  and  that 
the  transient  shear  stress  parameter  Af1 1  Fig.  5.13  is  larger  than 
that  given  in  Ref.  [4]  when  the  10-point  computational  molecule  with 
H  =  T  =  .01  is  used  in  the  solution  of  the  unsteady  momentum  equation. 
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CURVES  FOR  CHAPTERS  III  AND  IV 
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Fig.  3.4  Transient  Profiles  of  the 
Stream  Function  Increment  Af 
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Fig.  3.9  Comparison  of  Closed  Form  and  Recurrence 
Solutions  for  the  Enthalpy  Increment  Ag  for  g  =  2 


Fig.  3.10  Comparison  of  Empirical  Steady  jState  Enthalpy 
Distribution  with  the  Exact  Solution  for  Me  |*  iq  and  Pr  =  0.72 


Fig.  3.12  Satisfaction  of  the  Unsteady  Energy  Equation 
by  the  Empirical  Enthalpy  Distribution  for  an 
Adiabatic  Wall  with  MP  =  10,  Pr  =  0.72 
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Fig.  5.2  Enthalpy  Jump  for  an  Adiabatic  Wall 
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Fig.  5.5  Transient  Heat  Transfer  Parameter 
from  Assumption  g(n »0)  =  0  with  M  =  10 
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Fig.  5.9  Transient  Heat  Transfer 
Parameter  from  Empirical  Distribution 
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Fig.  5.11  Steady  State  Weak-Interaction 
Induced  Pressure  for  M  =  10,  x  =  4.5 
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APPENDIX  3 
FORTRAN  PROGRAM 


*•  IN  COL.  6:  THIS  LINE  IS  A  PART  OF  PREVIOUS  STATEMENT 


C  CALCULATION  OF  DELTA  F  DISTRIBUTION  USING  DOUBLE  PRLCI 

*  SION 

C  AND  USING  THE  TRANSFORMATION  X  I  =  T AU/ (  I f TAU ) 

EXTERNAL  FCT,OUTP 

DI MENS  I  ON  PRM T (5 ) . AUX (  16*3 )  »Y<  3)  * DERY( 3 )  ,F0{ 501), FOPRI 

*  M<  501 )  » 

IF0  2PR [ <  501  )  ,R ( 501  ) ,CC ( 50  1  )  . SI <  50 1  )  ,52(501  )  . S3( 50  l  ) 

REAL *8  NUM , PRM T , AUX, Y  , DERY , FO ♦ F OPR  I M , F 0 2PR I  , S 1  , S2, S3, S 

*  , R  »  CC , TOP  * 

1  DO T , EP S , H , T , DENOM , A , 8 , C , D , E , X  I 
COMMON  / BLOCK 1 /FO/ 8LOCK2/FO PR I M/ OL  OCK3/F02PR I/8LUCXA/N 
-X  N/8LOCK5/H 
EPS-0 .0 1  DO 
NN-50 1 
NNN=NN- l 
MM  =  5  1 
VMM -MM— 1 
H=  0 . 0  IDO 
T  —  0 .0100 

C  GENERATION  OF  8LASIUS  SOLUTION  USING  SSP  HAMMINGS  PRED 

*  ICTOR 

C  CORRECTOR  METHOD 

PRMT (  1)  —  0.0D0 
PRMT (2 )=5.0D0 
PRMT ( 3 ) -H 
PRMT ( 4 ) =0 . 5D-6 
ND I M=3 

DERY ( l ) = .33300 
DERY ( 2  )  -  .333D0 
DERY ( 3 )= .334D0 
Y ( 1 ) =0 • ODO 
Y ( 2 )  =  0 .0  00 
Y ( 3 ) — 0 . 4696 00 ODO 

CALL  DHPCG  (  PRMT  ,  Y  »  DERY  »  NO  I  M  *  IHLF  »FC  T  ,  GU  TP  ,  AUX  ) 

C  NUMERICAL  SOLUTION  OF  THE  MOMENTUM  TO  (  3  •  ^  )  UoINC,  RECU 

*  R PENCE 

C  RELATIONSHIP  (3.11) 

C  APPLYING  BOUNDARY  CONDITION  (2.42) 

S  1  (  l ) =0. ODO 
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c 


c 

c 

c 


c 


S2(1)-0.0D0 
S  3 (  1  )  =  0  .  ODO 

APPLYING  BOUNDARY  CONDITION  (2.42) 
on  l 0  N=1 ,NNN 

S  l  (  N  + 1  )  -  3  1  ( N  )  +  H* ( EP  S* (  1 • - F 0 PR  I M ( N )  )  ) 

10  CUNTINUE 

WR I TE( 6 .20) 

HH-0  .0 

DU  30  N  = 1.496*5 

WRI TE( 6 , 50)  HH, S 1 ( N) , S  1  <  N+  1  )  ,  S  1  ( N+2 )  *  SI  { NT  3 )  »  S  l  (  NT4  ) 
HH-HH+0 . 05 
30  CONTINUE 
HH— 5.0 

WR I T  E ( 6 , 50 )  HH , S 1 ( NN ) 

WRITING  DELTA  F  VALUES  ON  TAPE 
WRITE ( 4 , ERR-40 )  SI 

EVALUATION  OF  FIRST  TIME  LINE  USING  7-POINT  CGMPUTATIQ 
*  NAL  MOLECULE 
OF  FIGURE  3.1 
R ( 1 ) =0.0 DO 
R  (  2  )  =  0 . 2  50  0 
CC(  1  )  =  0 .0D0 
CC( 2 ) -0 .ODO 
R l=R (  1  ) 

R2=R (2  ) 

C 1 =CC (  1  ) 

C  2- CC ( 2 ) 

CALL  FI RST  {NN,H»T,Rl  , R2 . Cl  ♦ C2, FO ,  FOPR I M  »  F02PR  I  »  S 1 *S2) 


X  I  =T 

T AU-SNGL ( X  1/ (  1  .-X  I  )  ) 

WRITE (6, 60)  XI  *  T AU 
HH—  0 . 0 

DO  7  0  N  =  1  *  49 6*5 

WRI TE ( 6  *50 )  HH,S2<  N )  *S2<  N+ 1  )  »S2 (N  +  2 )  , S2< N  +  3 )  , S2 ( N  +  4 ) 
HH— HH+O .05 
70  CONTINUE 
HIT-5 .0 

WR I TE ( 6  *  50 )  HH»S2(NN) 

WRITING  DELTA  F  VALUES  ON  TAPE 
WRI TE( 4  ,ERR=40 )  S2 
DO  30  M=3*MM 


X I=X I +  T 

TAU  =  SNGL ( XI /(  1  .-XI  )  ) 

DO  90  N-3.NNN 

A-4.  O'AT  +  4 . 0*H*T*FO  (  N)  +6 .0*H*H*  (1  .O-XI  )*(X  I*F0PR  IM(N  )-  1 
*  .  O  +  XI  ) 

B=- 1 2 . 0*T-8 • 0 *H* T*FO ( N ) +4 .0*H*H*H*F 02PR I  (N  )  *<  f-3 . O'* X  I* 


*  (  1 . 0  -  X  I  )  ) 

C=12.0^T+4.0  X=H  *  T*F  0  (  N  )  +6  «  0  *  (  1 

*  OPRIM(M) ) ) 

D— —  4  •  0  T 

F  =  8  ♦  0  *  (  1  .O-XI  ) I  *  0  —  X  I  Y (  l 

*  2  (  N  +  1  >  ) 


O-XI  )  *H*H* (  1  • 0— X  I  * (  l  • OTF 


0  TF OPRIi-K  N)  )  )  Y  (  S2  (  N -  1  )  -S 


»  •  r 

. 

r 
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1  +  <  4 . 0*X  l  *<  1  •  0  — X  I  )  *H*H*H*F02PR I  { N)  ) * ( S 1 ( N >  —4 . 0* S2 ( N )  )  +< 
*  2  •  0* ( 

21 • O-XI ) *( 1 • O-X I ) *H*H )*(Sl(Nfl)-Sl(N-l))+(2.0*XI*{l. O-X 
*  I  ) 

3*H*H*F0PRIM { N) ) * ( S 1 ( N- l ) -S 1 (N+ t ) ) 

R( N ) ——A/ ( B  TC*R ( N— 1  )+D*R(N-l  )*R(N-2)  ) 

NUM=E- (C*CC { N- 1 > TO*R (N-2 ) *CC { N- 1 ) +  D*CC ( N-2 ) ) 
DENGM=3+C*R( N-l ) +0*R<N-1 ) *R (N-2) 

9  0  C  C ( N )  =  N  U M/D E NO M 

C  USING  THE  STARTING  POINT  (3.15) 

THP-CC (NN-1 )-CC<NN-2 )/(4 ,0-R(NN~2) ) 

BOT  =  3 . 0/ ( 4 . 0-R(NN-2)  )-R(NN-l ) 

S3 ( NN ) =  T  OP/BOT 
I  I  — NN— 2 
DO  100  1=1,11 

N=NN— I 

S  3 ( N ) =R { N ) *  S3 ( NT  1  ) TCC ( N) 

100  CONTINUE 

WRITE (6, 60)  XI,TAU 
HH=  0  ,0 

DO  110  N= 1,496,5 

WRI TE ( 6 , 50 )  HH , S3 ( N )  , S3 (NT  1  )  »S3 (N  +  2 )  ,S3 ( N  +  3 )  , S  3( NT 4 ) 
HH=HH NO  .  05 
1 l 0  CONT INUE 
HH=5 • 0 

WR I TE ( 6 , 50 )  HH , S  3 ( NN ) 

C  WRITING  DELTA  F  VALUES  ON  TAPE 

WRI TE< 4,ERR=40)  S3 
DO  120  N  =  2 , N  N 
S  1  (  N ) =  S  2  (  N) 

S  2  {  N  )  =  S  3  (  N  ) 

120  CONTINUE 
SO  CONTINUE 

20  FORM  AT {  IH1  ,  ’DELTA  F  DISTRIBUTION  FOR  XI=0.00  IE  FOR  T 

❖  AU  =  0 .00*//) 

50  FORMAT ( 1 H  , F t 0 . 2 , 5F2 2 . 8 ) 

60  FORM  AT (  l HI  »* CELT A  F  DISTRIBUTION  TOR  xi  =  *F6.3f*  IE  TO 

*  R  TAU  =*F6«3 
1//) 

STOP 

40  W R I TE( 6, 1 30 ) 

130  FORMAT ( *  ERROR  WHILE  WRITING  ON  TAPE') 

STOP 

END 


. 

■ 
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SUBROUTINE  FCHX.Y.DERY) 
DOUBLE  PRECISION  DERY(3),Y<3) 
DERY  (  1  )—Y(2.) 

DERY ( 2) =Y( 3) 

DERY ( 3)=-Y ( 1 ) *Y( 3) 

RETURN 

END 
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SUDROUT INE  OUTP (X,Y, DERY , 1 HLF , ND I M , PRMT ) 

REALMS  Y(3)  iFO (501  )  , FOPRIM (501  >  ♦ F02PR  I  (  501  )  , H 
COMMON  /BLOCK  1 /F0/BLGCK2/F0PR I M/BLOCK3/F02PRI/DLOCK4/N 
*  N/ BLOCK 5/ H 
C=0 . 0 
X  0=  0  .  0 
AH=SNGL ( H ) 

DO  10  J  —  1  »NN 
l  =  J 

X J=X0+C*AH 

I F ( ABS ( X  J—  X )  • LT • 1  «E  — 3 )  GO  TO  2  0 
1 0  C=C+1 .0 
RETURN 

20  F  0 (  I  ) — Y (  1 ) 

FOPR  IM  (  I  )  =  Y  (  2  ) 

F02PR I ( I )=Y( 3) 

RETURN 

END 


■ 
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SU3RQUT INC  FIRST  (NN,H»T»R1  ,R2(C1  ,C2trO, FOPR IM *  F  0  PP R I  ♦ S 

*  1  »  S2 ) 

C  CALCULATION  OF  FIRST  TIME  LINE  USING  RECURRENCE  SGLUTI 

*  ON  WITH 

C  7-POINT  MOLECULE  OF  FIGURE  3.1 

01  MENS  ION  FO ( 50 1  )  *  F 0 P R I M (501), F02PR I (  50 1  )  *  R { 5 0 l )  , CC ( 50 

*  1  )  .SI  ( 501  )  , 

IS 2 ( 50  1  ) 

REAL*8  NUM  .  DENOM  .  FO*  FOPR  IM  »  F  0  2  P  ft  I  ,  R  .  CC  ,  S  1  »  S2  ,  A  *  B  ,C  *  O  »  E 

*  ,  X  I  .  H  ,  T  ,  R  1  , 

IR2.C1  »C2 .TOP  *  BOT 
NNN=NN-1 

R(  1  )  -Rl 

R( 2 )=P2 

CC( 1 )~C1 

CC ( 2)=C2 
X  I=T 

DU  10  N=3.NNN 

A=T+H*T*F0<  N  (  1  .  0-X  I  )  *  (  1  •  0 -X  I  *  (  1  •  0  TFOPR  IM(  N  )  )  ) 

0=-3 .  0*T-2 .  0*H*T  *F0< N ) +  H*H *H*T *F 0 2PR I  ( N ) -2 • Q*X I  * (  1 . 0-X 

*  I ) *H*H*H* 

1F02PRI (N) 

C=3 . 0 *T +H*T*F0 ( N ) +H*H*(  1  .0-X I  ) * (  l  •  0-X I *(  1  •  0  TF  OPR  I M ( N ) ) 

*  ) 

D  T 

E=H*H*  (1  •  0  -  X  I  )  *  (  1  ■»  0  — X  I  *(  1 . 0  +  FCPR  IM(  N  )  )  )  *(  S  1  (N-l  )  -SI  (  N+ 
*  1  )  >- 

12.0*XI*<  1 . 0-X I  ) *F02PR I ( N ) *H *H* H* S 1  {  N) 

R (N ) =— A/ ( B +C*R(N-1 ) +D*R(N- 1 ) *R ( N-2 ) ) 

NUM=E- (C*CC (N-l ) +D*R ( N-2 } *CC (N- 1 ) +D*CC( N-2 ) ) 
DENOM=D+C*R< N-l ) +  D*R  <N- 1 )*R(N-2) 

CC( N ) -NUM/OENOM 
10  CONTINUE 

TOP=  CC  (  NN-  1  )  -CC  (  NN-2  )  /  (  4 .0-R  (N.N-2)  ) 

BOT  =  3 . 0/ ( 4 .0-R( NN-2 )  ) -R ( NN- 1  ) 

S  2 ( NN ) -TOP/OQT 
I  I ~ NN-2 
DO  20  1=1,11 

N=NN— I 

5  2 ( N )=R(N)*S2(N+1  ) + CC ( N ) 

20  CONTINUE 
RETURN 
END 


r  • 
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appendix  c 

FORTRAN  PROGRAM 


'**  IN  COL.  6:  THIS  LINE  IS  A  PART  Oh  PREVIOUS  STATEMENT 


C  SOLUTION  OR  UNSTEADY  STATE  ENERGY  EQUATION  ASSUMING  GD 

*  OT=0  AT  TAU= 

C  CASE  OF  MACH  NUMBER  =  10*  PR  AND  TL  NUMBER  =  0.72  AND  ISO 

*  THERMAL  WALL 

EXTERNAL  AFCT , FC T , DFC T , OUTP 

REALMS  MACH* GAMMA , PR  *  H*  T , NUM, DENOM *  U  *  GS  *  GSPR I M ♦ GS2PR I  * 

*  P.FPRIM, 

1 F2PR IM, FOOT ,  X  I  , A » 8 , C , D. R * CC *F0 *FC2PR.I  , GO  »F OPR  I M.F3PRI 

*  M »  CC C  »  BB  *  RR , 

2Y »  AUX.PRMT .A A.DERY ,FDPRI M* GW  »  GODOT  , MEO *  EPS , G1 
D I  MENS  ION  U { 3 . 50 1  )  , GS ( 50 1  )  , GSPR I M ( 50 1  ) , GS2PR I  ( SO  1 )  , R{ 5 

*  0 1 ) *CC { 50 1 ) , 

IFO (501)  *  F02PR 1(501)* G0( 50 I)  *0(501  )  ,FPRI  M  (  50  I )  .F2PRIM (5 

*  01), FOOT (501 

2 )  , FOPR IM ( 50  1  )  *  PRMT ( 5)  , F3PR I M (50 1  )  ,CCC ( 2  *  2 )  , 88 ( 2 , 2 ) , Y ( 2 

*  ) , AUX( 20 , 2) , 

3RR{ 2)  , AA( 2 ,2 )  *  OERY ( 2 )  , FDPRIM ( 50 1  )  , GODOT  (  50  1  )  ,  G  1  (  501  ) 
COMMON/BLOCK l /F 
COM MON /BLOCK 2/ F2  P  R I M 

C OM MON / SL OCK3/FOPR I M/ BLOCK4/FO 2PR  I 

C  O M MON/ B  L  OC  K  5/G 1 

MACH=1 0 .0D0 

PR=0. 7200 

GAMMA= 1 . 4D0 

GW=0 .000 

NN-50 1 

NNN=NN-t 

M  M  =  51 

.MMM  — MM  —  1 

H=0 .0100 

T=O.OIDO 

EPS  =0 .0100 

C  FVALUATING  MACH  NUMBER  CORRESPONDING  TO  INITIAL  STEADY 

*  STATE 

ME0=DSQRT { M ACN*M ACH/ ( 1 • +MAC H#MACH #EPS * ( GAMMA— 1 • ) ) ) 

C.  READING  THE  81.  AS  I  US  VALUES  STORED  ON  TAPE 

READ(  3  ,ERR-250 )  FO 
RE  A  O ( 3  *  ERR=  2  5  0 )  FOPRIM 
R  E  A  O ( 3 , ERR-250 *FND= 10)  F02PR I 


■ 
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10  CONTINUE 

C  EVALUATING  STEADY  STATE  G 2  ENTHALPY  DISTRIBUTION  FROM 

*  SOLUTION  (3. 

C  USING  SSP  SIMPSON’S  RULE 

C  ALL  F  I  N AL I  PR ,  1  ♦ GW ,  MACH , GAMMA ,  F 0 2 PR  I  ,  GS  ) 

W R  I  T  E { 6  *  2  0 )  PR  *  M A  CM 

20  FORMAT ( 1H1 * ’STEADY  STATE  G  DISTRIBUTION  FUR  PR=  »F6.2, 

*  *  AND  MACH 

1  NO •  *  F6 • 2  »  *  USING  SUBROUTINE  FINAL’//) 

HH= 0*0 

DO  1  N=1 ,496,5 

WRITE!  6,30)  HH  ,  GS{  N  )  ,  GS(  N+  1  )  ,  GS  (  iNf  2  )  ,  GS  (  N  +  3  )  ,  GS(  N+4  ) 
HH=HH+0.O5 
1  CONTINUE 
HH-5 .0 

WRITE (6, 30)  H  H , G  S { N  N ) 

C  WRITING  STEADY  STATE  VALUES  ON  TAPE 

WRITE <4 , ERR =40)  GS 

C  EVALUATION  OF  FIRST  AND  SECOND  STEADY  STATE  DERIVATIVE 

*  5  USING  SSP 

CALL  DDET  5 ( H , GS , GSPR I M , NN ,  IER) 

CALL  DDET5 ( H , GSPR I M , GS2PR I , NN ♦ TER ) 

C  READING  MOMENTUM  EQUATION  SOLUTION  FROM  TAPE 

DO  50  N  = 1  , NN 

READ ( 3 , ERR=2  50 )  FIN)  , FPR I M ( N )  , F2PR I M ( N )  , FDO  TIN) 

50  CONTINUE 

FDPR  IM (  l  )=  < -3 . *F  DOT { 3 ) +4 . *FDO  T I  2 ) -FOOT (  1  )  ) / I  2 • *H ) 

DO  60  N=2,NNN 

FOPR IM ( N )  =  I  FOOT ( N+ l )-FDOT IN- I  )  >/(2.*H) 

60  CONTINUE 

FDPR IM I NN )= I  3. ❖FOOT I NN ) -4 . *FDQT I NN- 1  ) +  FDOT I NN-2 )  )  /  I  2  •  ❖ 

*  H ) 

C  USING  DIFFERENTIAL  EQUATION  TO  OBTAIN  THIRD  DERIVATIVE 

C  USING  VECTOR  SPACES  FOPRIM  AND  F02PRI 

DO  70  N— 1 , NN 

F3PR  I  M  I  N  )  —2  •  ❖FDPR  I  M  IN)  -FIN)  *F2PR  IM  I  N  ) 

FOPR I M ( N ) =F2PR I M ( N ) 

F02PR I  I N ) =F3PR IM IN ) 

70  CONTINUE 

C  EVALUATION  OF  G1  DISTRIBUTION  USING  SSP  DLBVP 

ND IM-2 

8B I l , 1 ) =1 . 0D0 
BBI 1 ,2) -0.0D0 
RBI  2,1  ) =0 . 0D0 
OL3I  2,2)  =0. 0D0 
CCC I  1  ,  1  ) =0  •  0D0 
CCC ( 1 ,2) =0 .ODO 
CCC I  2, 1  ) =1  .00  0 
CCC I  2  *  2 )-0 .ODO 
RRI l ) =  GW 
RR ( 2 )  =  1  .ODO 
DERY I  1  J-0.5D0 


' 


■ 


- 


•  ■ 4 
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DER Y ( 2 ) =0.500 
PRMT  ( 1  ) =0 .000 
P  R M T  (2 ) =5.0  00 
PRMT ( 3 ) =H 
PRMT ( 4 ) =0 . 50-4 

CALL  DLBVP( PRMT , 80  ,  CCC  »  RR , Y , DFRY , ND I M ,  IHLF , AF  CT  * FCT , DF 
*  CT iOUTP »AUX* 

1  AA  ) 

WR I TF ( 6*80)  IHLF 

80  FORMAT (  I  HO  *  1  ERROR  PARAMETER  IHLF  IS  *  I  5// ) 

WRITE (6, 90) 

90  FORMAT!  1H1 ,  *  INI T TAL  <31  DISTRIBUTION  USING  SUBROUTINE  D 
4  LBVPV/) 

HH  =  0 . 0 

OO  100  N= 1 » 496, 5 

WRI TE ( 6 , 30 )  HH, G L { N ) »  G 1  <  N  +  1  )  , G1  { N  +  2 )  , G1 ( N  +  3)  , G 1 ( N+4) 

HH=HH+0 . 05 
100  CONTINUE 
Hi  1  =  5.0 

WRITE! 6.30)  HH»G1 ( NN) 

C  USING  SUBROUTINE  FINAL  TO  CALCULATE  INITIAL  STEADY  STA 


C 


c 


*  TF  ENTHALPY 
D  ISTR I BUT  I ON 

CALL  F I NAL ( PR ,  I  , GW , ME 0 , GAMMA , F  0  2PR I  *  GO ) 

WRI TE(6. 1 10) 

110  FORMAT (  1  HI  ENTHALPY  JUMP  FROM  ASSUMPTION  GDOT  =  0’//) 
DO  120  N — 1 • NN 
GODOT ( N ) =G 1 ( N ) -GO ( N ) 

120  CONTINUE 


HH=0 .0 

DO  130  N= 1,496, 5 

WRITE (6, 30 ) HH, GODOT ( N) , GODOT! N FI ) , GODOT (N+2) , GODOT ( N+3 
A  )  ,  GO D  OT  {  N  +4  ) 

HH—HH+0 .05 
130  CONTINUE 


HH=5.0 

WRITE  (6, 30)  HH,  GODCJT  {  NN  ) 

EVALUATING  INCREMENTAL  ENTHALPY  DISTRIBUTION  USING  GOD 
4  OT  VECTOR  SP 
X  1=0 .0 
T  AU  =  0 . 0 


WR  I  T  E  (  o  ,  1  4  0  )  X  I  ,  T A U 

140  FORM  AT (  1  HI ,  ' DELTA  G  DISTRIBUTION  FOR 


X  I  =  '  F  6  •  3  ,  ' 


IE  F 


*  OR  TAU=  * 

IF6 . 3//) 

DO  150  N=1,NN 

U( 1 » N) =GS ( N ) -G1 ( N) 

U  (  2  »  N  )  =  U  (  1  »  N  ) 

GODOT (N)=GS(N)-U( 1 ,N)-G1(N) 
150  CONTINUE 


HH=  0 . 0 

DO  1 60  N=t  »  496 ♦ 5 


•  . 


- 


* 


- 


' 
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WR I TE( 6 , 30 ) HH, GODOT ( N)  ,  GODOT ( N  +  1  )  , GODOT { N  +  2 )  ,  GODOT ( NT 3 

*  )  *  GODO  T ( N+4) 

HH= HH  +  O .05 

160  CONTINUE 
HH=5 • 0 

WRITE! 6,30)  HH, GODOT (NN) 

X  I=T 

T  A  U  =  S  N  GL ( X I/(  1  .-XI  )  ) 

WRITE! 6, 140)  XI , TAU 
HH=  0 . 0 

DO  l 70  N= 1 , 496,5 

WRIT  E( 6 *30) HH, GODOT ( N ) , GODOT ! N  +  1  )  , GODOT ( N  +  2 )  , GODOT <  N  + 3 
=5=  ),  GODOT  (N  +  4) 

HH=  HH  +0 .05 
170  CONTINUE 
HH=5. 0 

WR I T  £ ( 6 , 3  0  )  HH , G  ODO T ( NN ) 

3  0  FORMAT! 1H  »  F l 0 . 2 , 5F  20.8) 

C  WRITING  INCREMENTAL  ENTHALPY  DISTRIBUTION  ON  TAPE 

WRITE! 4  *  ERR  =  40 )  (U(l.N)  ,N=1,NN) 

WRITE! 4  , ER R= 4 0 )  ( U ( 2 , N )  , N~ 1  , NN ) 

C  SOLVING  THE  UNSTEADY  STATE  ENERGY  EQUATION  FOR  TAU  >  0 

DO  130  N=1,NN 

READ!  3 »  ERR= 2  50 )  F{N)*FPRIM(N)  ,F2PR I M (N)  , FOOT ( N) 

130  CONTINUE 

R ( 1 )=0 .ODO 
CC<  l  ) =  0 • OD  0 
DO  190  M~3 »  MM 
XI-X I +  T 

TAU=SNGL( X I/(  1  .-XI  )  ) 
on  200  N=1 ,NN 

READ ( 3 , ERR— 250 )  F ( N )  ,FPRIM(N)  ,F2PR IM(N)  , F  DOT (N) 

200  CONTINUE 

DO  210  N=2  * NNN 

A- T * !  I  . /PR+H* ( F( N ) /2 .-XI  * (  l  .-XI) *FDQT ( N )  )  ) 

Q— -3 . H* (  1  . -X  I  * (  1  .  +FPR I M ( M )  )  ) * !  1  .-XI  ) - 2 • * T/PR 
C  =  T* (  1  • /PR— H* (F ( N) /2 . -X  I* ( 1  .-XI) *F DOT ( N )  )  ) 

D=— H^H * ( 1 .-XI *( l • +FPRIM!N) ) ) * ( 1 .-XI ) *(4.*U(2»N)-U! 1 *N) 

*  ) +M*H*T* ( 

1 GS2 PR  I  ! N ) /PR+ ! F ( N ) -2 • *X I  * ( 1  .-X  I  ) ❖FOOT ( N )  ) ❖GSPR  I  M ( N ) + ( G 

*  A  M  M  A  —  1  .  )  * 

2MACH*MACH*F2PR IM ( N ) *F2PR IM (N ) ) 

DENOM— B+C^RIN— 1 ) 

R ( N ) =— A/DENOM 
NUM=D-C*CC(N-1 ) 

C  C  (  N  )  .=NUM/DENOM 
210  CONTINUE 

U( 3 ,NN) -0. ODO 

GODO  T ( NN ) =GS ! NN ) -U ( 3 , NN ) -G 1  ! NN ) 

C  MARCHING  DOWNWARD  USING  THE  RECURRENCE  RELATION 

DO  220  I = 1 , NNN 
N=MN-I 


- 
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U ( 3.N)=R(N)*U(3,N+1) +  CC ( N ) 

GO  DOT { N ) =GS ( N ) -U ( 3  »  N ) -G 1 ( N ) 

22  0  CONTINUE 

WRITE{6» 140)  XI* TAU 
HH-0 .0 

DO  230  N  =  1  , 496,5 

WR IT  E( 6 , 3  0 ) HH, GODOT ( N)  *  GODOT (N  +  l  )  *  GODOT ( N+2 )  , GODOT! N+3 
*  ) , GODOT ( N  +  4 ) 


HH=HHfO .05 

230 

CO NT  INUE 

HH— 5.0 

WRITE! 6,30)  HH , 

G 00 OT ! NN) 

V; R  I  T  F  (  4, ERR=40  ) 
00  240  N—  1  ,  NN 

U( l ,N)-U(2,N) 

U  (  2  ,  N  )  -U  C  3  ,  N  ) 

!  U  !  3  , 

II 

Z 

21 

NN  ) 

240 

CONT I NUE 

190 

CONTI NUE 

ENDFILE  4 

STOP 

250 

WRI TE ( 6, 260 ) 

260 

FORMAT ( *  ERROR 

STOP 

WHILE 

READING 

FROM  TAPE' ) 

40 

WRITE! 6, 2 70 ) 

270 

FORMAT! »  ERROR 

W  Hi  I  LF 

’WRIT  I  NG 

ON  TAPE*) 

STOP 

END 

* 
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SUBROUTINE  AFCT  (  X  ■»  A  A  ) 

REAL *8  AA( 2* 2 ) , F( 501 > *PR 

COMMON/BLOCK 1 /F 

PR-0 . 7200 

A  A (  1,1  )=0.0D0 

A  A (  1  , 2 ) =1  .00  0 

AA( 2, 1  )  —  0.000 

NN=50 l 

H=.01 

C=0 .0 

X0  =  0 .0 

DO  10  T  =  1  ,  N N 
N=I 

XI-X0+C*H 

IF ( ABS( XI-X)  .LE.O  .005)  GO  TO  20 
C=C+1 . 

10  CONTINUE 


WRITE ( 6* 30)  XI , X 

30  FORMAT ( 1  HO,  • AFCT  COULD  NOT  MATCH 
STOP 

20  AA( 2 ,2) =-PR*F (N) 

RETURN 

END 


XI 


AND  X='2F15.5) 


SUBROUTINE  FCT(XfFF) 

REAL*8  PR  ,  GAMMA , MACH, F2PR I M ( 50  1 )  , FF { 2 ) 

COMMON/DLOCK2/F2PRIM 

PR=0 • 7200 

GAMMA- 1,400 

MACH-1  0,000 

FF (  1  )  =0. 000 

N  N  =  501 

H=,01 

C=0 . 0 

XO-O  .0 

DO  10  I =1 , NN 
N—  I 

X I=X0+C*H 

IF ( ABS( X I-X ) ,LF . 0 .005 )  GO  TO  20 
C  =  C  +  1  . 


10  CONTINUE 

WRITE! 6, 30)  XI,X 

30  FORMAT (  1  HO*  • FCT  COULD  NOT  MATCH  XI  AND  X= * 2F 1 5 . 5 ) 
STOP 

20  FF( 2  )=-PR*( GAMMA -1  .  ) *MACH*MACH*F2PR IM{ N ) *F2PR IM (N ) 
RETURN 
END 
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SUBROUTINE  DFCT(X.DF) 

RF.AL*8  OF  (  2  )  ,FOPR  IM  (  50  1  )  ,  F 0  2 PR  I  (  SO  1  )  ,PR  ,  GAMMA  ,  MACH 

COMMON/DL OCK3/FOPR I M/BLGCK4/F02PR I 

GAMMA=! . 400 

PR=0 . 72D0 

MACH- l 0 .000 

OF {  1  )=0 .000 

NN=501 

H=  .0  1 

C^O  .0 

xo-o.o 

DO  10  I  =  1 , NN 
N=  I 

X I=XO+C*H 

1F( ABS( XI-X)  .LE .0  .005 )  GO  TO  20 

c-c+i . 

10  CONTINUE 

WRITE! 6,30)  Xf,X 

3  0  FORMAT (  1  HO ,* DFCT  COULD  NOT  MATCH  XT  ANO  X=I2F15.5) 
STOP 

2  0  Or ( 2  )  =  -2 •♦PR*! GAMMA- 1 .  ) *MACH*M ACH*FOPR I M{ N ) *F02PR l  (  N ) 
RFOURN 
END 


. 


SU3RC-UT  INt  OUTP(X,V,  DERV  ,  I  HLF  *  ND  I  M  ,  PRMT  ) 
RF  AL*3  Y ( 2)  , G1 { 501  ) 

COM MGN/BLOC K  5/G 1 

N  N  =  50 1 

H-  .  0  1 

C=0 . 0 

X0=0 .0 

DO  10  1=1 iNN 

N=  I 

X  I =XC  +  C*H 

IF ( ABS{  XI— X)  •  L  T ♦ 0 • 5E  —  3 )  GO  TO  20 
C  =  C  +  1  . 

10  CONTINUE 
20  G1(N  )  =  Y(  1  ) 

RETURN 

END 
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c 

c 


c 

c 


c 


c 


c 


c 

c 


SUL3ROUT  I  NE  FI  NAL  (  PR  ,  I  ,  GW  ,  MACH,  GAMMA  ,  F2PR  IM  ,  GS  ) 

1=0  FOR  ZERO  WALL  HEAT  TRANSFER  SOLUTION 

1=1  FOR  NON-ZERO  WALL  HEAT  TRANSFER  SOLUTION 

DIMENSION  R  <  50 1  )  ♦ AF ( 50 1  )  , AF  1  ( 50  1  ) , F2PR I M ( 50 1  )  , S { 50  1  )  , G 

*  S{ 50  1  ) 

REAL  *8  PR  •  H,  MACH,  GAMMA  ,  R  ,  AF  ,  AF  1  ,  F2PR  I M  ,  S  ,  G S  ,  E  ,  G  VV 
NN  =  50 1 
H=  0 . 0  IDO 

INTEGRATION  BY  SSP  SIMPSON'S  RULE 

EVALUATION  OF  AF ( N )  USING  AF 1 { N )  VECTOR  SPACE  TEMPORAR 

*  ILY 

F  =  2 . -PR 
DO  10  N= l , NN 

I  F  (  F 2 P R  I M  (  N  )  *  L  E  •  0  •  0  )  GO  TO  20 
AF 1  ( N )  =  ( F2 PR  I M ( N )  ) * *E 
GO  TO  10 
2  0  AFI  (N)=0  .0 
10  CONTINUE 

CALL  DQSF ( H , AF 1 , AF , NN) 

EVALUATION  OF  AF 1 ( N ) 

DO  30  N=1,NN 

I F ( F2PR  I  M <N  )  .LE . 0 . 0 )  GO  TO  40 
AFI  ( N ) =  (  ( F2  PR  I M ( N  )  ) * *PR ) 4 AF { N ) 

GO  TO  30 
4  0  A  F 1  ( N  >  =  0  .  0 
30  CONTINUE 

EVALUATION  OF  R ( N  )  USING  AF ( N )  VECTOR  SPACE  TEMPORARIL 

*  Y 

CALL  DOSE (H, API , AF,NN) 

DO  50  N=1 ,NN 

R ( N ) =2 . * PR* ( AF { NN ) - AF ( N ) ) 

50  CONTINUE 

•  IF (  I  .EQ  .  1  )  GO  TO  6  0 

EVALUATION  OF  STEADY  STATE  G  DISTRIBUTION  FOR  ZERO  WAL 
4  L  HEAT  TRANS 
DO  7  0  N  = 1  *  N N 

GS ( N ) = 1 , + ,5* ( GAMMA- 1 . ) ❖MACH*M ACH*R ( N ) 

70  CONTINUE 
RETURN 
60  CONTINUE 

EVALUATION  OF  STEADY  STATE  G  DISTRIBUTION  FOR  NON-ZERO 

*  WALL  HEAT  T 

EVALUATION  OF  S(N)  USING (N)  AND  AF ( N )  VECTOR  SPACES 
DO  80  N= 1 ,NN 

IF(F2PRIM(N)  .LE.0.0)  GO  TO  0  0 
API  ( N ) = ( F2P  R  I  M ( N )  ) * *PR 
GO  TO  80 
90  AF 1 ( N ) =0 . 0 
80  CONTINUE 

CALL  DOSE  (H*  AFI  •  AF.NN.) 

DO  100  N=1  »  NN 

S  (  N  )  =  AF  (  NN  )  -  AF  {  N  ) 


( 
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GS  (  N  )  =  l  .  +.5*( GAMMA-1  .  )  *M  ACH*MACH*R  (  N  )-( 1 .+.5*<  GAMMA-1  . 
*  ) *MACH*MACH 

1  * R (  1  )  —  GW ) *S ( N ) /S ( 1  ) 

100  CONTINUE 
RETURN 
END 


' 
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APPENDIX  C 
FORTRAN  PROGRAM 


» *  i 


IN  COL.  6:  THIS  L  I  NT  IS  A  PART  OF  PREVIOUS  STATEMENT 


C  CALCULATION  OF  UNSTEADY  STATE  HEAT  TRANSFER  PARAMETER 

C  AND  TRANSIENT  INDUCED  PRESSURE  WITH  WALL  HEAT  TRANSFER 

C  CASE  OF  MACH  NO. =1 0 , PRA NDTL  NO. =0 . 72 , CHI =  4 . 5  AND  GW  =  0 

*  (COLD  WALL) 

REA  L* 3  NU , GW , GPR I M  » G , GS , U , FO PR  I M , FO  2PR I , H  , GS  PR  I M , U  PR  I M 

*  tNUSt 

IGAWSt  MACH, PR , X , GAMMA. DNU , KE I, Y, FP RIM, I, J, 12, 01 ,0J, Yl , Z 

*  1  »  T  ,  P  »  A  »  X I  , 

2Z,II,J1,DU,DJ1,I21,  NU  I  ,GPRIM1  ,  DNU  1 ,  P  1  ,  P  2  ,  E  ,  Yt  *  GO  PR  I  M  , 
-  N  E  t  ON  E  ,  IE, 

3  J  E  »  D I  E  »  C J  E , R  E  C , GA  WS 1 »  NU  S 1 , G  $  PR  I  1,  I  F?,MEO 

DIMENSION  NU( 51 ) , G ( 50 1 ) , GS ( 50 1 ) , U ( 50 1 ) , FOPR I M { 50 1 ) , GPR 

*  I M (51). 

1F02PR l ( 501 ) ,NU1 ( 51  )  ,  CN'J(51  ), 1(51), J(51), 01(51), DJ (5 

*  1 ) ,Y( 501 ) , 

2FPR IM( 501 ) , Yl { 501 ) , Z1 ( 501 ) , P( 51 ) , Z ( 5Q1  )  ,  I  1 ( 51  )  ,  J1  ( 51  )  , 

*  011(51), 

3DJ 1 ( 5  1  > ,DNU1 ( 51 ) ,P1( 51 ) ,P2( 51 )  ,E( 5  01)  , YE( 501 )  ,NE(  51 ) ,0 

*  NE(51), 

4  I E ( 5 1 ) ,JE(51)»DIE(51) ,DJE(51) 

MACH=1C.CDQ 
PR  =  0. 7200 
GW=0. ODO 
CM  I  =  4 .500 
GAMMA =1. 400 
T=0.01D0 
H  =  0 .0  IDO 
MM  =  5  1 
N  N  =  5  0  1 

REODSQRT  (  PR  ) 

C  EVALUATION  OF  ADIABATIC  WALL  ENTHALPY  RATIO 

G  A  W  S  =  1  .  +  R  E  C  * 0 . 5  *  (  G  A  A  M  A  -  1 .  )  *  M  AC  H*  M  A  C  H 
G A W S  1  =  1  . 4-0 . 5*1  GAMMA-1 .  )  *M ACH*M ACH 
A=DSORT (2. DO) 

C  READING  BLASIUS  VALUES  E ROM  TAPE 

R  E  A  D  (  3  ,  EP  R  =  1  0  )  (  X  ,  N=  1  ,  N  N ) 

R  F  A  D  (  3  ,  t  P  R  •=  1  0  )  (FCPRIM(N)  ,  N  =  1  ,  N  N ) 

READ( 3 , EP  R  =  1 0 , E  N  0  =  2  0 )  (F 02  PR  I (N)  , N  =  1 , N N ) 

20  CONTINUE 


ir 


m 

' 


■ 
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C 

C 

c 


c 


c 


0 

c 

c 


c 


c 


READING  STEADY  STATE  GS  VALUES  WITH  WALL  HEAT  TRANSFER 
R  E AD ( 4  »  ER  R  =  1 0 )  GS 

EVALUATION  OF  STEADY  STATE  INTEGRAL  12 
INTEGRATING  USING  SSP  LIBRARY  SIMPSON'S  RULE 
DO  30  T  =  1  f  N N 
Y1(N)=GS(N)-F0PRIM(N) 

Y( N  )=  1 . +0 . 5*( GAMMA-1 .  ) *MACH*MACH* FOPR  IN ( N ) * ( 1 . -FOPP I M ( 

*  N))+(GW-1.)* 


1  (  1. -FOPP ININ)  )  -  F  0  P  R  I  M  (  N  ) 

USING  VECTOR  SPACE  E  TEMPORARILY  FOR.  EMPIRICAL  ENTHAL 
*  PY  D I S  Tp I BU  T 


E ( N )  =  1 .+( GW-1 .  ) *( 1 . -FOPP  IM(N  )**3)  +  ( PR**.32  3  )*< 1 . *. 5*<G 

*  ATM A- 1 • ) *MAC 

1H*MACH*CS0RT ( PR  )  -  GW ) *FGPR I M ( N ) * { l . -FOPR  IM { N  )**2 )- . 5*< G 

*  A  M  M  A  - 1  .  )  *  P  R * 

2  (  FOPP  IP  (N  )  *  *  2  )  *  (  l.-FOPFIM(N)  )  *M  ACH*  .VAC  H- FOPR  I M  (  N  ) 

USING  VECTOR  SPACE  G  TEMPORARILY 
G(N) =  Y< N) +  E0PR I M ( N ) 


30 


❖ 


CONTINUE 

CALL  DCSF(H,YltZl ,TN) 

I  2  -  Z  1  {  N  N  ) 

CALL  DCSF(H,Y,Z,NN) 

I  2  1  =  Z  (  N  N  ) 

CALL  CCSF(H.E. Z» NN) 

I  E  2  =  Z  (  N  N ) 

EVALUATING  STEADY  STATE  SLOPE  AT 
ORWARD  DIFFE 

r'  S P R  I  M  =  (-11.  *GS  (  1  )  +1  8  .  *GS  (  2  )  -9  .  * 
EVALUATION  OF  STFADY  STATE  HEAT 


PLATE  SURFACE  USING  F 

GS ( 3  )  +  2. *GS( 4 ) ) /( 6.*H ) 
TRANSFER  PARAMETER 


N  L)  S  =  G  $  P  R I M / (GAtoS-GW) 

GSPR  I  1  =  (-1 1 .*G ( 1 ) +18. *G ( 2 )-9.*G( 3 ) +2.*C( 4 )  ) /( 6 . *H ) 

NUS  1  =  GS  PR  II / ( GANS1 -GW ) 

EVALUATION  OF  INTEGRALS  I  AND  J  BY  SIMPSON’S  RULE 
00  40  P=1,MM 
DO  50  N  =  1 , NN 

R  E  AD ( 3 , ER  R  =  1 0  )  X , F  PR  I M ( T  )  , X  ,  X 
50  CONTINUE 

READING  DFLTA  G  (=U)  VALUES  FROM  T A P E ( W I T H  WALL  HEAT  T 


*  RAN  SEEP) 

P  E  A  D ( 4  ,  E R R - 1 0 )  ( U ( N ) . N= 1 » NN ) 


DO  60  N  =  1 »  NN 

G(N)=l.+.5*( GAMMA-1. ) *MACH*M ACH*F P R  I M  (  N  ) 

*  +(GW-1.)*<1. 

1  -  F  P  R  I  M  (  N  )  ) 


(  l.-FPR  I M ( N ) ) 


Y1 ( N ) =GS( N ) -U ( M J-FPRIM (N) 

E  (  N  )  =  1 .  +  (  GW-  1 .  )  *  (  l.-FPPIM(N) 
*  MMA-1.)*MACH 
1*m  ACH*C  SOP  T  i  PR  ) - G W  )  *'FPR  I  M  (  N  ) 


* ❖ 3 )  +■  ( PR** .  33  3  )*( 1 

*( 1 .  -  FPR I M ( N ) *  * 2  )  - 


*  A  -  1  .  )  *  P  R  *  (  F  P 

2R  I  m  (  M  )  i?  if  2  )  *(  1  •  -FPR  IM  (  N  )  )  *M A CM* MACH 
USING  Y  VECTOR  FOR  PR  =  1  SOLUTION 


+■  .  5  *  (  G  A 
5*( GAMM 


Y(N)=G(N)-FPRIM(N) 


' 
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V  L-  !  N  )  =  E  (  N  )  -  F  P  P  I  N  (  N  ) 


USING  FORWARD  DIFFER 


H) 


GPRIM(M)=GSPRIM-UPRI M 


H) 

H) 


C  EVALUATION  OF  FEAT  TRANSFER  PARAMETER 


M 1 1  (  M ) =  G PR  I M 1 / ( GA W S 1 -G W ) 


DNU1  (  M  )  =  N  U  1  (  M  )  -NUS1 
N  U  {  M )  =  G  PR  I  N  (  M  )  /  (  G  A  W  S -G  W ) 

DNU ( M ) - N  U ( M ) - N U  S 
N  E  (  M  )  =  G  E  P  R  I M  /  (  G  A  W  S  -  G  W  ) 

D.ME  (  M)=NE(  MJ-NUS 
CALL  C-  C  S  F  (  H ,  Y  1  »  Z  t  N  N  ) 

I ( M )- Z  (NN) 

CALL  OCSF(H,Y,Z,NN) 

1 1  (  M  )  =  Z  (  N  N  ) 

CALL  DCSF(H,YE,Z,NN) 

I E ( M ) - Z ( NN) 

DO  7 C  N  -  1  *  N N 

Y  1  {  N  )  =  G  S  {  N  )  - 1 J  (  N  )  - 1  • 

Y  (  N  )  =  C  (  N  )  -  l . 

Y  E  (  N  )  =  E  (  N1  )  - 1  . 

70  CONTINUE 

CALL  DCSF(H,Y1,Z1,NN) 

J(M)=Z1 (NN) 

CALL  0CSF(H,Y»Z1 »NN) 

J 1  ( M ) =  Z I ( NN) 

CALL  CCSF(H,YE»Z1 » N  N  ) 

J  E  (  M  )  =  Z  1  (  N  N ) 

40  CONTINUE 

WR I T  E ( 6 , 80 )  PR , M ACH , GW , NUS 

80  FORMAT (  1H1 ,' STEADY  STATE  HEAT  TRANSFER  PARAMETER  FOR  P 

*  R='F6.2,'  MA 

1CH  NO*  '  F  6  •  2  *  '  GW  'F6.2,'  IS  'FIG.  5//) 

WRITE (6 ,90) 

90  FORMAT  (  1 H  , 'UNSTEADY  HEAT  TRANSFER  PARAMETER  NIJ'//) 

WR I T E ( 6 , 1 00 )  ( NU ( M ) , M= 1 » MM ) 

100  FORMAK  1H  ,5F26.8) 

WRITE (6, 110) 

110  FORMAT (1H  , »  I  NCR EMENTAL  HEAT  TRANSFER  PARAMETER'//) 

WR I T  E ( 6 ,1 00 )  ( DNU ( M ) , M  =  !  , MM ) 

WRITE (6, 120)  GW 

120  FORMAT ( 1HC. 'UNSTEADY  GPP  I M  FOP  GW  =  ’F6.2//) 

WRITE! 6 , 1 3 0 )  GW , G  S PR  I M 

130  FORMAT!  LHO, 'STEADY  STATE  GSPRIM  FOR  GW=  'F6.2,  '  IS  'FI 

*  0.4//) 

WRITE (6, 140)  PR 

140  FORM AT ( 1H1 » ' HEAT  TRANSFER  PARAMETER  FROM  EMPIRICAL  DIS 

*  TR  I  BUT  ION  FC 


■ 

' 
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1R  PR  =  ’ F6 .2/ /  ) 

WR  I  T  F  (  6 , 1 00  )  (  NE  (  M  )  ,  M=  1  ,  MM  ) 

WRITE (6, 110) 

W  R I  T  F  (  fc  ,  1 C  0  )  ( DME ( M  )  , M= 1  , MM ) 

WRITE! 6*150) 

150  FORMAT (  1HL , 1  FEAT  TRANSFER  PARAMETER  FRGV  CLOSED  F 0 R M  S 

❖  CLUTICN  FOR 
iPR=  1'//) 

WRITE (6, 160) 

160  FORMAT! 1H  ,» UNSTEADY  HEAT  TRANSFER  PARAMETER  NU1  FROM 

❖  CLOSED  FORM 
1  SOLUTION'//) 

WRI T E ( 6  *  1 00 )  (  MJ  1 ( M )  , M= 1 , MM ) 

WRITE! 6, 110) 

W  R I TF ( 6 , 1  GO )  ( DNU 1 ! M )  , M  =  1 , MM ) 

WRITE (6, 170)  12,121 

170  FORMAT! 1H1 STEADY  STATE  INTEGRALS  I?,  121  ARE  »2r15.3/ 

❖  / ) 

WRITE (6 ,180) 

160  FORMAT! IHO, *  I TEGRALS  I  *  I  1  ,  J  ,  J  1  ARE  '//) 

WR I T  E ( 6 , 1 00 )  (  I  ( N  )  ,  M= 1  ,  MM ) 

WRITE (6, 100)  !  II  (M) , M- 1 , MM) 

WRITE (6 ,100)  ( J(V)  ,  M  = 1 , MM ) 

WR I T  F ( 6  ,  1 00 )  !  J  1  ( M  )  , M= 1  , MM ) 

C  EVALUATION  OF  INTEGRAL  DERIVATIVES  01  AND  DJ 

CALL  CD  El  5(T,  I  ,01  ,f  4,  IFF  ) 

C  ALL  DC E T  5  (  T  ,  J  ,  0  J  ,  f-  M  ,  I  E F  1  ) 

CALL  CDET5 ( T ,  I  1 , D 1 1 , MM , I E R ) 

CALL  CPET  5 ( T , J I , D J 1 , MM ,  I  f R  1  ) 

CALL  DPET5(T, IE, DIE, MM,  I  ERE) 

CALL  DDET  5 ( T , JE • D JE  »  N  M , I E  RE ) 

WRITE (6,190) 

190  FORMAT! IHO,' DER  I  VAT  IVES  P I  ,  D 1 1 , D J , 0 J 1  A*E  •//) 

WR I T  E ! 6 , ICO )  (PI ( M  )  , M  = 1 , MM ) 

WR I T  E ( 6 , 1 00 )  ( D I  I ! M ) , M= 1 . MM ) 

WR I TE ( 6 , 1 OC )  ( 0  J ( M )  , M  = 1 , MM ) 

WR I T E ( 6 , ICO )  ( D J 1 ( M ) , M= 1  *  MM ) 

C  EVALUATION  CF  TRANSIENT  INCUCEQ  PRESSURE 

K  1=0 
K  =0 

200  CONTINUE 
X  I =C. OCO 
DO  2  10  M  = 1 , M M 

P  (  M  )  =  !  GAMMA  /  (  MACh*MACH  )  )  *  (  I  2/A-A*  (  1 .  -  X  I  )  *  (  1  .—XI  )  *DJ  (  '-1 ) 

❖  -(l./A)*(I(M 

1  )-2,-XIv(  1  .-XI  )  *  D  I  ( M )  )  ) 

p i ( M )  =  ( GAMMA/ ( MACH^MACH )  ) * ( A* ( 1 . -X  I ) * (  1  .-XT  ) *0J ! M )  M  1  • 

❖  /  A  )  *  <  I  (  M  )  -  2  . 

l*X  I  *  (  1  .  -X  I  )  *D  KM))) 

X  I =X I +T 
210  CGNTINUF 

WRITE (6, 220)  PR 


' 


. 

* 
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220  FORMAT (  1H1 ,*  INDUCED  PRESSURE  FUNCTION  (P/PE) /CHI  FOR  P 

*  R= ' F6 • 2/ / ) 

WRI TE(6 ,100)  (  P(  M)  ,  fv  =  1  ,  MM) 

WRITE (6,230)  PR 

230  FORMAT (  1H0 INDUCED  PRESSURE  FUNCTION  NOT  REFERENCED  T 

*  0  FINAL  S  T  E A 

ID Y  STATE,  FOR  PR=’F6.?//) 

WR  I  T  E  (  6 , 1  00  )  (  P  1  m  ,  1  ,  MM  ) 

C  EVALUATION  OF  INDUCED  PRESSURE  FOR  CHI =  4.5  AND  MACH  NO 

*  *  =  10 

C  (USING  VECTOR  SPACE  J) 

DO  2 AO  M=1  , MM 
J (M) =CF  I*P(M!) 

P  2  (  M  )  =C  H I  *  P 1  (  M  ) 

240  CONTINUE 

WRITE (6,250)  CHI, MACH 

250  FORMAT* 1H1 , « INDUCED  PRESSURE  FOR  CHI  PARAMETER ' F6 .2, * 

❖  AND  NACH  NO. 

1  *  F  6 . 2  /  /  ) 

UR  I  TE ( 6  ,  ICC  )  ( J{ M  )  ,  M=1 , MM ) 

WRITE (6,260)  CHI , MACH 

260  FORMAT (  ihO  ,*  INDUCED  PRESSURE  NOT  REFERENCED  TO  FINAL  S 

❖  TEADY  STATE, 

IFOR  CHI  PARAMETER* F6.2  ,  •  AND  MACH  N0.'F6*2//) 

WR I T E ( 6 , 1 00 )  ( P 2 ( M )  ,  M= I , MM ) 

I  F  (  i<  1  « C-  E  .  I  )  STOP 
IFfK.GE. 1 )  GO  TO  270 
DO  2  80  8  =  1,  MM 
I  (  M  )  =  I  E  (  M  ) 

D I  ( M ) =D  I  F { M ) 

DJ  (  M  )  =  D  J  E  (  M  ) 

280  CONTINUE 
K  =  K  +  1 
GO  TC  200 
2  70  K  L  =  K  1  +  3. 

PR=  1 . 

DO  29  0  M  = 1  ,  M M 
I  (  M  )  =  1 1  (  M  ) 

D I  (  M ) =D 1 1 ( M ) 

DJ ( M)=DJ1 ( M) 

290  CONTINUE 
12=121 
GO  TO  200 
10  WRITE (6,300) 

300  FORMAT (  *  ERROR  WHILE  R  E  A  D I N G  FROM  TAPE') 

STOP 

END 


. 

. 


